Strange density behavior with tr. octahedral box and constant pressure

From: Michael Jakusch <mjakusch_at_pharma.anbi.ethz.ch>
Date: Fri 7 Sep 2001 15:24:53 +0200

Dear AMBER users,

I'm somehow confused by the output I get when running constant
pressure MD with sander (v6) using a truncated octahedral box. I hope,
somebody will be able to comment on what I did wrong ...

I prepared my system using leap:

  -----snip-----
  ....
  .. ( buildt unit "tk_oct" ) ..
  ....
> solvateOct tk_oct WATBOX216 9.0 0.8
  Scaling up box by a factor of 1.261587 to meet diagonal cut criterion
    Solute vdw bounding box: 79.771 70.800 66.919
    Total bounding box for atom centers: 102.480 93.509 89.628
    Solvent unit box: 18.774 18.774 18.774
  The number of boxes: x= 6 y= 5 z= 5
  Adding box at: x=0 y=0 z=0
  Center of solvent box is: 46.935873, 37.548698, 37.548698
  ....
  .. ( 149 more boxes )
  ....
    Total vdw box size: 105.242 96.666 92.522 angstroms.
    Volume: 470627.482 A^3 (oct)
    Total mass 248603.900 amu, Density 0.877 g/cc
    Added 9538 residues.
  ....
  -----snap-----

So far, everything looks reasonable for me since the calculated volume
is exactly the half of the volume of the cuboid with the given box
size, and the density is in the usual range.

However, when I started the NPT calculation (after performing several
preceding equilibration steps), I get this output:

 -----snip-----
 ....
 NSTEP = 0 TIME(PS) = 6.500 TEMP(K) = 300.82 PRESS = -186.93
 Etot = -81934.0770 EKtot = 25131.9596 EPtot = -107066.0366
 BOND = 1686.7582 ANGLE = 4707.5129 DIHED = 6431.5103
 1-4 NB = 2332.3271 1-4 EEL = 23119.0222 VDWAALS = 7403.2870
 EELEC = -152746.4543 EHBOND = 0.0000 CONSTRAINT = 0.0000
 EKCMT = 9241.8855 VIRIAL = 13040.7453 VOLUME = 941254.9631
                                                Density = 0.4386
 Ewald error estimate: 0.9434E-04
 ....
 -----snap-----

Here the value of the volume has doubled while the density is at half
the original value and is gradually increased to >~ 1 in the course of
the simulation!

I first suspected, that sander was erroneously treating the system as
a standard rectangular box, but from the header of the MDOUT it seems
to recognize everything correctly (I attached the full MDOUT-header
below).

The box information on the last lines of the PRMTOP and the previous
RESTRT files also seems ok for me:

 prmtop:
  0.00000000E+00 1.05242172E+02 9.66657070E+01 9.25220070E+01

 restrt:
 105.2421720 96.6657070 92.5220070 90.0000000 90.0000000 90.0000000


Thank you in advance for any suggestions

Michael


BTW: sander seems to _always_ read the box info from the file INPCRD,
regardless of the value of NTX



-- 
Dr. Michael Jakusch
ETH - Swiss Federal Institute of Technology
Departement of Applied Biosciences
Winterthurerstrasse 190
CH-8057 Zürich
Switzerland
Phone: +41.1.635 60 71
Fax:   +41.1.635 68 84
email: michael.jakusch_at_pharma.anbi.ethz.ch
*****************************************************************
Header of MDOUT
*****************************************************************
          -------------------------------------------------------
          Amber 6  SANDER                   Scripps/UCSF 1999
          -------------------------------------------------------
|         Thu Sep  6 18:20:26 2001
File Assignments:
|MDIN : _sandin_                                                              
|MDOUT: tko_npt1.out                                                          
|INPCR: tko_nvt2.crd                                                          
|PARM : tk_apo_oct.prmtop                                                     
|RESTR: tko_npt1.crd                                                          
|REFC : refc                                                                  
|MDVEL: mdvel                                                                 
|MDEN : tko_npt1.mden                                                         
|MDCRD: tko_npt1.mdcrd                                                        
|MDINF: mdinfo                                                                
 Here is the input file:
#L1                                                                            
#L2                                                                            
#L2                                                                            
 &cntrl                                                                        
  imin=0,                                                                      
  ntx=7                                                                        
  nstlim=50000,dt=.002		                                                       
  ntpr=10,ntwe=10,ntwx=100	                                                    
                                                                               
  ntc=2,ntf=2,tol=0.00001	                                                     
  ntb=2, cut=9, scee=1.2                                                       
  ntt=1, tautp=.2		                                                            
  ntp=1, taup=.05		                                                            
				                                                                           
				                                                                           
  temp0=300			                                                                 
  ibelly=0                                                                     
 &end                                                                          
-------------------------------------------------------------------------------
#L1                                                                             
   1.  RESOURCE   USE: 
 getting box info from bottom of parm
 getting new box info from bottom of inpcrd
| peek_ewald_inpcrd: Box info found
   EWALD SPECIFIC INPUT:
 -------------------------------------------------
 NO EWALD INPUT FOUND: USING DEFAULTS
 -------------------------------------------------
     Largest sphere to fit in unit cell has radius =    46.261
     Calculating ew_coeff from dsum_tol,cutoff
     Box X =  105.242   Box Y =   96.666   Box Z =   92.522
     Alpha =   90.000   Beta =   90.000   Gamma =   90.000
     NFFT1 =  108       NFFT2 =   96       NFFT3 =   96
     Cutoff=    9.000   Tol   =0.100E-04
     Ewald Coefficient =  0.30768
     Interpolation order =    4
 NATOM  =   39522 NTYPES =      17 NBONH =   34484 MBONA  =    5150
 NTHETH =   11434 MTHETA =    7042 NPHIH =   21356 MPHIA  =   13262
 NHPARM =       0 NPARM  =       0 NNB   =   94936 NRES   =   10488
 NBONA  =    5150 NTHETA =    7042 NPHIA =   13262 NUMBND =      43
 NUMANG =      87 NPTRA  =      39 NATYP =      31 NPHB   =       1
 IFBOX  =       2 NMXRS  =      24 IFCAP =       0
   EWALD MEMORY USE:
|    Total heap storage needed        =       1641
|    Adjacent nonbond minimum mask    =      94936
|    Max number of pointers           =         25
|    List build maxmask               =     189872
|    Maximage  =      57820
   EWALD LOCMEM POINTER OFFSETS
|      Real memory needed by PME        =       1641
|      Size of EEDTABLE                 =      20768
|      Real memory needed by EEDTABLE   =      83072
|      Integer memory needed by ADJ     =     189872
|      Integer memory used by local nonb=    2635590
|      Real memory used by local nonb   =     647724
|    MAX NONBOND PAIRS =   40000000
|     Memory Use     Allocated         Used
|     Real            20000000      2274558
|     Hollerith        4800000       247622
|     Integer         80000000      4342895
|     Max Nonbonded Pairs:40000000
     BOX TYPE: TRUNCATED OCTAHEDRON
   2.  CONTROL  DATA  FOR  THE  RUN
                                                                                
     TIMLIM=  999999.   IREST =    0       IBELLY=    0
     IMIN  =    0
     IPOL  =    0
     NTX   =    7       NTXO  =    1
     IG    =    71277   TEMPI =     0.00   HEAT  =    0.000
     NTB   =    2       BOXX  =  105.242
     BOXY  =   96.666   BOXZ  =   92.522
     NTT   =    1       TEMP0 =  300.000
     DTEMP =    0.000   TAUTP =    0.200
     VLIMIT=    0.000
     NTP   =    1       PRES0 =    1.000   COMP  =   44.600
     TAUP  =    0.050   NPSCAL=    0
     NTCM  =    0       NSCM  = 9999999
     NSTLIM=50000       NTU   =    1
     T     =    0.000   DT    =   0.00200
     NTC   =    2       TOL   =   0.00001  JFASTW =    0
     NTF   =    2       NSNB  =   25
     CUT   =    9.000   SCNB  =    2.000
     SCEE  =    1.200   DIELC =    1.000
     NTPR  =      10    NTWR  =      50    NTWX  =     100
     NTWV  =       0    NTWE  =      10    IOUTFM=       0
     NTWPRT=       0    NTWPR0=       0    NTAVE=       0
     NTR   =    0       NTRX  =    1
     TAUR  =   0.00000     NMROPT=    0       PENCUT=   0.10000
     IVCAP =    0       MATCAP=    0       FCAP  =    1.500
   OTHER DATA:
     IFCAP =    0       NATCAP=    0       CUTCAP=    0.000
     XCAP  =    0.000   YCAP  =    0.000   ZCAP  =    0.000
     VRAND=    0
     NATOM =   39522  NRES =  10488
     Water definition for fast triangulated model:
     Resname = WAT ; Oxygen_name = O   ; Hyd1_name = H1  ; Hyd2_name = H2  
   3.  ATOMIC COORDINATES AND VELOCITIES
     Largest sphere to fit in unit cell has radius =    46.261
 NEW EWALD BOX PARAMETERS from inpcrd file:
     A     = 105.24217  B    =  96.66571  C     =  92.52201
     ALPHA =  90.00000  BETA =  90.00000  GAMMA =  90.00000
                                                                                
 begin time read from input coords =     6.500 ps
 Number of triangulated 3-point waters found:     9814
     Sum of charges from parm topology file =  -0.00000028
     Forcing neutrality...
 ---------------------------------------------------
 APPROXIMATING switch and d/dx switch using CUBIC SPLINE INTERPOLATION
 using   5000.0 points per unit in tabled values
 TESTING RELATIVE ERROR over r ranging from 0.0 to cutoff
| CHECK switch(x): max rel err =   0.3302E-14   at   2.581700
| CHECK d/dx switch(x): max rel err =   0.8064E-11   at   2.761360
 ---------------------------------------------------
     Total number of mask terms =      84636
     Total number of mask terms =     169272
|    Total Ewald setup time =   0.15999985
 ------------------------------------------------------------------------------
Received on Fri Sep 07 2001 - 06:24:53 PDT
Custom Search