CHARMM Development Project Forums User Discussion & Questions Molecular Dynamics incorrect density when using VV2 to run drude

Hi,
I was trying to reproduce the simulation of a paper by Anisimov et al. in 2007 (Polarizable Empirical Force Field for the Primary and Secondary Alcohol Series Based on the Classical Drude Model)
I set up a 45.8x45.8x45.8 (A^3) box which consists of 128 molecules. (this initial density is about 10 times smaller than the experimental value, and hence the volume should decrease in the NPT simulation to reach a correct density)
I ran a NPT simulation for 250 ps at 298K and 1 atm with VV2 integrator, but the equilibrium volume is way too high (the final box length is 160 A)(it should decrease not increase).
The box length didn't increase nonstop though, it reached 160 then fluctuated around it for the rest of the simulation.
I didn't get any error message from CHARMM, but this result seems suspicious ( although the final report average temperature and pressure is fine). I am new to CHARMM so I am not sure what part of setting I got wrong.
I already tried a bunch of different setting such as removing shake constraint, remove the cutoff in crystal build command, remove the image command, etc. Nothing worked.
If I increased the target pressure, the equilibrium volume will decrease which suggests the barostat is working. I've also run a similar simulation with the FF in lammps (but I still wanted to run it in CHARMM), the equilibrium density is very close the experimental value( as is the case in the paper).
Any help would be really appreciated

Here's my inp file:
set toppar toppar
stream @toppar/toppar_defs_2018.str
! read the psf and coordinate file
read psf card name meoh.psf
read coor card name meoh.crd


coor sdrude
coor shake

! shake constraint
SHAKE bonh param nofast -
select .not. type D* end -
select .not. type D* end

! SETUP CRYSTAL (DEFINE, BUILD),
set boxval = 45.8
crystal define cubic @boxval @boxval @boxval 90. 90. 90.
crystal build noper 0 cutoff 12


! Set up images -- center the protein by segment and the solvent by residue
image byres xcen 0.0 ycen 0.0 zcen 0.0 sele all end


! set up nonbond parameters
nbond nbxmod 5 VDW vatom vswitch e14fac 1.0 cutnb 14.0 cutim 14 ctofnb 12.0 ctonnb 10.0 -
inbfrq -1 imfrq -1 wmin 1.0 -
ELEC atom cdiel EWALD PMEWALD KAPPA 0.34 FFTX 32 FFTY 32 FFTZ 32 ORDER 6
energy


tpcontrol nther 2 nstep 20 -
ther 1 tau 0.1 tref 298.15 select .not. type d* end -
ther 2 tau 0.005 tref 1 select type d* end -
baro btau 0.1 pref 1.00

set NSTEPS = 20000
set printfreq = 1000
OPEN WRITE CARD UNIT 62 NAME test.nose ! info for temperaature
open unit 35 write form name test.rst ! restart file
open unit 32 write unfo name test.trj ! trajectory
open unit 33 write form name test.kunit ! Temperature
open unit 34 write unfo name test.vel ! velocity

DYNA vv2 start -
nstep @NSTEPS time 0.001 -
iunread -1 -
iunwrite 35 iuncrd 32 iunvel 34 kunit 33 -
nsavc @printfreq nsavv @printfreq nprint @printfreq isvfrq @printfreq -
iprfrq @NSTEPS ntrfrq @NSTEPS - ! what's used for swm4 water
firstt 298.15 -
iasvel 1 iseed 314159 -
iuno 62 nsnos @printfreq

output result:
* * * AVERAGES FOR THE LAST 250000 STEPS
AVER DYN: Step Time TOTEner TOTKe ENERgy TEMPerature
AVER PROP: GRMS HFCTote HFCKe EHFCor VIRKe
AVER INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
AVER EXTERN: VDWaals ELEC HBONds ASP USER
AVER IMAGES: IMNBvdw IMELec IMHBnd RXNField EXTElec
AVER EWALD: EWKSum EWSElf EWEXcl EWQCor EWUTil
AVER PRESS: VIRE VIRI PRESSE PRESSI VOLUme
---------- --------- --------- --------- --------- ---------
AVER> 250000 250.00000 3228.26570 684.18595 1956.40578 298.14828
AVER PROP> 15.75889 2640.59173 0.00000 587.67397 644.22682
AVER INTERN> 204.05459 292.06749 9.81547 38.87825 0.00000
AVER EXTERN> -2.00684 350.69952 0.00000 0.00000 0.00000
AVER IMAGES> -0.12272 -6.81850 0.00000 0.00000 0.00000
AVER EWALD> 45.90840 -97862.98641 98886.91653 0.00000 0.00000
AVER PRESS> -34.2523 -395.2322 0.5524 1.0296 4206262.9857
---------- --------- --------- --------- --------- ---------
Lattice Parameters> Averages for the last 250000 steps:
Crystal Parameters : Crystal Type = CUBI
AVER A = 160.67117 B = 160.67117 C = 160.67117
AVER Alpha = 90.00000 Beta = 90.00000 Gamma = 90.00000
AVER PIXX = 0.97 PIYY = 1.09 PIZZ = 1.04
AVER PIXY = -0.06 PIXZ = -0.05 PIYZ = -0.05
AVER Gradient Norm = 2.70061



The toppar_defs_2018.str file is from charmm/toppar/drude, I just copy all the files to another folder named toppar
the meoh.psf and meoh.crd are obtained from another separated script:
set toppar toppar
stream @toppar/toppar_defs_2018.str

open unit 91 card read name initial.pdb
read sequ pdb unit 91
set residue MEOH

generate @residue first none last none setup warn drude dmass 0.4

rewind unit 91
read coor pdb unit 91

ic param
coor sdrude
coor shake
coor print


open write unit 10 card name meoh.crd
write coor card unit 10

open write unit 10 card name meoh.psf
write psf card unit 10

open write unit 10 card name meoh.pdb
write coor pdb unit 10

stop

where initial.pdb file is just a 45.8x45.8x45.8 box with 128 meoh molecules (without drude particle and lone pairs)
the crd file itself looks normal (no undefined coordinates), here's part of it:
1280 EXT
1 1 MEOH C1 17.1060000000 5.5580000000 26.9880000000 MEOH 1 0.0000000000
2 1 MEOH DC1 17.1060000000 5.5580000000 26.9880000000 MEOH 1 0.0000000000
3 1 MEOH O1 16.0640000000 6.2290000000 26.2880000000 MEOH 1 0.0000000000
4 1 MEOH DO1 16.0640000000 6.2290000000 26.2880000000 MEOH 1 0.0000000000
5 1 MEOH HO1 15.4600000000 5.5380000000 26.0950000000 MEOH 1 0.0000000000
6 1 MEOH H1A 16.9130000000 4.6950000000 27.6730000000 MEOH 1 0.0000000000
7 1 MEOH H1B 17.8380000000 6.2450000000 27.5060000000 MEOH 1 0.0000000000
8 1 MEOH H1C 17.6460000000 5.0320000000 26.2260000000 MEOH 1 0.0000000000
9 1 MEOH LP1A 15.8241442420 6.3464524504 26.5142170136 MEOH 1 0.0000000000
10 1 MEOH LP1B 16.1342348244 6.2343440269 25.9451611008 MEOH 1 0.0000000000
11 2 MEOH C1 39.1960000000 40.6790000000 26.7270000000 MEOH 2 0.0000000000
12 2 MEOH DC1 39.1960000000 40.6790000000 26.7270000000 MEOH 2 0.0000000000
.......



Thanks in advance!

One possibility is that your 125 MeOH molecules are all close together in the middle of your 45.8 cubic box, and that therefore no, or too few, images are generated. The output should tell you if this is the case; you would expect 26 images. Try with smaller box, it should be sufficient with a couple of extra Å on each side to avoid initial clashes.

Hi,
thank you for your reply, this seems to be the problem.
I only have 9 images instead of 26

: updating the image atom lists and remapping
Transformation Atoms Groups Residues Min-Distance
1 C001 has 30 3 3 2.73
2 C002 has 340 34 34 0.00
3 C003 has 50 5 5 0.63
4 C004 has 50 5 5 0.00
5 C005 has 440 44 44 0.00
6 C006 has 90 9 9 2.50
7 C007 has 400 40 40 0.00
9 C009 has 230 23 23 0.00

Is this the information regarding images?
I tried a smaller box, but when I ran the simulation, the volume still increased to about 160^3 A^3.
I also tried another initial configuration where molecules are placed neatly with equal distances (previously I used the insert-molecules function in gromacs to generate my initial pdb file).
In addition, I tried using nonpbc to first minimize my structure.
But I still ended up with a way too large volume with any initial configurations I have.
I also tried changing the IMGFrq to 50 instead of -1 and set IXTfrq to 500, still didn't work.

After CYRSTAL BUILD the numnber of transformations (images) is printed:
The number of transformations generated = 26

Do you get 26 with your smaller box? This is essential.
I do not understand why you think you need such a large box. If there are problems with a box of the correct size (eg high vdW energies due to clashes with images) it is usually sufficinent to increase the box size with 1Å.

Building a much smaller box from scratch and first use nvt to relax it seems to solve the problem.
Thank you for pointing out the source of the problem.
So does the number of images always have to be 26?
How do I know if my images are proper for the different systems?
I read the documentation regarding image but I still don't quite get it.

For parallelpiped unit cells, i.e. those with six parallel faces, the minimum number is 26, which is one image per face, edge, or vertex. The unit cell is at the center of a 3x3x3 array, completely surrounded by 26 image cells. (Crystal unit cells with edges shorter than the nonbond cutoff require additional images to satisfy the nonbond list building.) The exceptions to 26 are typically the more spherical shapes, such as the rhombic dodecadedron and the truncated octahedron, which require fewer images to completely surround the unit cell, e.g. 12 for the dodecahedron.