AMBER Archive (2001)

Subject: Strange density behavior with tr. octahedral box and constant pressure

From: Michael Jakusch (mjakusch_at_pharma.anbi.ethz.ch)
Date: Fri Sep 07 2001 - 08:24:53 CDT


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 ------------------------------------------------------------------------------