Thread Like Summary
lin1209
Total Likes: 2
Original Post (Thread Starter)
incorrect density when using VV2 to run drude #37883 05/18/2020 1:43 AM
by lin1209
lin1209
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!
Liked Replies
Re: incorrect density when using VV2 to run drude #37884 May 18th a 05:57 PM
by lennart
lennart
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.
1 member likes this
Re: incorrect density when using VV2 to run drude #37886 May 19th a 07:06 AM
by lennart
lennart
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Å.
1 member likes this
Powered by UBB.threads™ PHP Forum Software 7.7.4
(Release build 20200307)
Responsive Width:

PHP: 5.6.33-0+deb8u1 Page Time: 0.368s Queries: 8 (0.365s) Memory: 0.8148 MB (Peak: 0.8449 MB) Data Comp: Off Server Time: 2020-09-29 08:09:34 UTC
Valid HTML 5 and Valid CSS