From: Batuhan Kav (bkavlist_at_gmail.com)
Date: Tue Mar 10 2020 - 05:08:14 CDT
Dear NAMD users,
I am facing problems with my pure POPC bilayer simulations using the
Charmm drude polarizable force field in NAMD. The main issue is that, as
soon as I use a time step larger than 0.1 fs, in the equilibration
phase, the simulations crash with the following error:
WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 1000
ERROR: Constraint failure in RATTLE algorithm for atom 4123!
ERROR: Constraint failure; simulation has become unstable.
I've been trying to simulate a POPC bilayer using the Charmm Drude
polarizable force field. The system consists of 36 POPC lipids in each
membrane leaflet. There are no ions involved. I solvated the system with
the SWM4-NDP water model.
I initially minimized and equilibriated the structure using Charmm36
force field. I have extracted the last frame of 100 ns simulation to use
with the Drude prepper in Charmm-gui.
I used the Drude prepper in Charmm-gui to obtain the coordinate files as
well as the NAMD input files. For the minimization, I have created a
constraint file with VMD in which I constrained the non-drude particles
using a scaling factor 10 (constraintScaling 10). With this input, I
have minimized the system for 2000 steps.
After the initial constrained minimization, I removed the constraints
from the system and re-minimized for another 2000 steps.
For the initial equilibration setup at 303.15 K, I used a flexible cell
while keeping the ratio constant. Furthermore, I used a time step of 0.1
fs while using the drudeDamping 20.0. I constrained bonds in the system
using "rigidBonds all" option. I also set the drudeTemp to 1. I
equilibriated the system for 100000 steps (10 ps)
Above equilibration step runs without any errors. However, I do observe
that within a few hundred steps the temperature of the drude particles
rises to 2.5 K (drudebond column in the log file) and although this
value fluctuates, it never goes down to 1 K, which I provided with the
drudeTemp option.
When I switch to the production run with 1.0 fs time step, within a
couple of hundred steps the temperature of the drude particles go beyond
10 K and the simulation eventually crashes with the above stated error.
I am also attaching the relevant parts of the input files I have used below.
Is there something I am doing wrong here? I am aware that the error
message is related to the improper minimization/equilibration, but I was
expecting that the system will equilibriate with 0.1 fs after 100.000
steps. However, even with a very small time step I don't seem to be
dealing with the Drude particles properly.
Any help would be greatly appreciated.
Thanks in advance,
Batuhan
1) Initial minimization with restraints
structure ../step2_drude.psf
coordinates ../step2_drude.pdb
# Force-Field Parameters
paraTypeCharmm on;
parameters ../toppar_drude/toppar_drude_master_protein_2013f.str;
parameters ../toppar_drude/toppar_drude_lipid_2017c.str;
parameters ../toppar_drude/toppar_drude_model_2013f.str;
parameters ../toppar_drude/toppar_drude_carbohydrate_2018a.str;
parameters ../toppar_drude/toppar_drude_nucleic_acid_2017c.str;
# These are specified by CHARMM
exclude scaled1-4
1-4scaling 1.0
switching on
vdwForceSwitching no;
mergeCrossterms on
# You have some freedom choosing the cutoff
cutoff 12.0;
switchdist 10.0;
pairlistdist 16.0;
stepspercycle 20;
pairlistsPerCycle 2;
# Integrator Parameters
timestep 0.1;
rigidBonds all;
nonbondedFreq 1;
fullElectFrequency 1;
# Use Drude polarizable model
drude on
drudeTemp 1
drudeHardwall on
drudeDamping 20.0
drudeBondLen 0.2
drudeBondConst 40000
drudeNbtholeCut 5.0
exec tr "\[:upper:\]" "\[:lower:\]" < ../step2_drude.str | sed -e "s/
=//g" > step2_drude.str
source step2_drude.str
cellBasisVector1 $a 0.0 0.0;
cellBasisVector2 0.0 $b 0.0;
cellBasisVector3 0.0 0.0 $c;
cellOrigin 0.0 0.0 $zcen;
wrapWater on;
wrapAll on;
wrapNearest off;
# PME (for full-system periodic electrostatics)
PME yes;
PMEInterpOrder 6;
PMEGridSpacing 1.0;
# Pressure and volume control
useGroupPressure yes;
useFlexibleCell yes;
useConstantRatio yes;
langevin on
langevinDamping 1.0
langevinTemp $temp
langevinHydrogen off
# constant pressure
langevinPiston on
langevinPistonTarget 1.01325;
langevinPistonPeriod 50.0;
langevinPistonDecay 25.0;
langevinPistonTemp $temp
constraints on
consref restraint.pdb
conskfile restraint.pdb
conskcol B
constraintScaling 10
minimize 2000
run 200
2) Initial equilibration with 0.1 fs
# Force-Field Parameters
paraTypeCharmm on;
# These are specified by CHARMM
exclude scaled1-4 # non-bonded exclusion policy to
use "none,1-2,1-3,1-4,or scaled1-4"
# 1-2: all atoms pairs that are
bonded are going to be ignored
# 1-3: 3 consecutively bonded
are excluded
# scaled1-4: include all the
1-3, and modified 1-4 interactions
# electrostatic scaled by
1-4scaling factor 1.0
# vdW special 1-4 parameters in
charmm parameter file.
1-4scaling 1.0
switching on
vdwForceSwitching no; # force-based switching of vdW
should not be used for Drude FF
mergeCrossterms on
# You have some freedom choosing the cutoff
cutoff 12.0; # may use smaller, maybe 10.,
with PME
switchdist 10.0; # cutoff - 2.
# switchdist - where you start
to switch
# cutoff - where you stop
accounting for nonbond interactions.
# correspondence in charmm:
# (cutnb,ctofnb,ctonnb =
pairlistdist,cutoff,switchdist)
LJcorrection yes
pairlistdist 16.0; # stores the all the pairs with
in the distance it should be larger
# than cutoff( + 2.)
stepspercycle 20; # 20 redo pairlists every ten steps
pairlistsPerCycle 2; # 2 is the default
# cycle represents the number of
steps between atom reassignments
# this means every 20/2=10 steps
the pairlist will be updated
# Integrator Parameters
timestep 0.1; # fs/step. If experiencing
instability, try smaller value, e.g., 0.75 or 0.5.
# If drudeHardwall is not
applicable, please set 0.5.
rigidBonds all; # Bound constraint all bonds
involving H are fixed in length
nonbondedFreq 1; # nonbonded forces every step
fullElectFrequency 1; # PME every step
# Use Drude polarizable model
drude on
drudeTemp 1
drudeHardwall on # Only available in latest NAMD
(newer than 2.09). Please comment out this line in old namd.
drudeDamping 20.0
drudeBondLen 0.25 # Please set 0.20 if
drudeHardwall is not applicable.
drudeBondConst 40000
drudeNbtholeCut 5.0
wrapWater on; # wrap water to central cell
wrapAll on; # wrap other molecules too
wrapNearest off; # use for non-rectangular cells
(wrap to the nearest image)
# PME (for full-system periodic electrostatics)
PME yes;
PMEInterpOrder 6; # interpolation order (spline
order 6 in charmm)
PMEGridSpacing 1.0; # maximum PME grid space / used
to calculate grid size
# Pressure and volume control
useGroupPressure yes; # use a hydrogen-group based
pseudo-molecular viral to calcualte pressure and
# has less fluctuation, is
needed for rigid bonds (rigidBonds/SHAKE)
useFlexibleCell yes; # yes for anisotropic system
like membrane
useConstantRatio yes
langevin on
langevinDamping 5
langevinTemp $temp
langevinHydrogen off
# constant pressure
langevinPiston on
langevinPistonTarget 1.01325
langevinPistonPeriod 50.0
langevinPistonDecay 25.0
langevinPistonTemp $temp
run 100000
3) Production with 1.0 fs
# These are specified by CHARMM
exclude scaled1-4 # non-bonded exclusion policy to
use "none,1-2,1-3,1-4,or scaled1-4"
# 1-2: all atoms pairs that are
bonded are going to be ignored
# 1-3: 3 consecutively bonded
are excluded
# scaled1-4: include all the
1-3, and modified 1-4 interactions
# electrostatic scaled by
1-4scaling factor 1.0
# vdW special 1-4 parameters in
charmm parameter file.
1-4scaling 1.0
switching on
vdwForceSwitching no; # force-based switching of vdW
should not be used for Drude FF
mergeCrossterms on
# You have some freedom choosing the cutoff
cutoff 12.0; # may use smaller, maybe 10.,
with PME
switchdist 10.0; # cutoff - 2.
# switchdist - where you start
to switch
# cutoff - where you stop
accounting for nonbond interactions.
# correspondence in charmm:
# (cutnb,ctofnb,ctonnb =
pairlistdist,cutoff,switchdist)
LJcorrection yes
pairlistdist 16.0; # stores the all the pairs with
in the distance it should be larger
# than cutoff( + 2.)
stepspercycle 20; # 20 redo pairlists every ten steps
pairlistsPerCycle 2; # 2 is the default
# cycle represents the number of
steps between atom reassignments
# this means every 20/2=10 steps
the pairlist will be updated
# Integrator Parameters
timestep 1.0; # fs/step. If experiencing
instability, try smaller value, e.g., 0.75 or 0.5.
# If drudeHardwall is not
applicable, please set 0.5.
rigidBonds all; # Bound constraint all bonds
involving H are fixed in length
nonbondedFreq 1; # nonbonded forces every step
fullElectFrequency 1; # PME every step
# Use Drude polarizable model
drude on
drudeTemp 1
drudeHardwall on # Only available in latest NAMD
(newer than 2.09). Please comment out this line in old namd.
drudeDamping 20.0
drudeBondLen 0.25 # Please set 0.20 if
drudeHardwall is not applicable.
drudeBondConst 40000
drudeNbtholeCut 5.0
wrapWater on; # wrap water to central cell
wrapAll on; # wrap other molecules too
wrapNearest off; # use for non-rectangular cells
(wrap to the nearest image)
# PME (for full-system periodic electrostatics)
PME yes;
PMEInterpOrder 6; # interpolation order (spline
order 6 in charmm)
PMEGridSpacing 1.0; # maximum PME grid space / used
to calculate grid size
# Pressure and volume control
useGroupPressure yes; # use a hydrogen-group based
pseudo-molecular viral to calcualte pressure and
# has less fluctuation, is
needed for rigid bonds (rigidBonds/SHAKE)
useFlexibleCell yes; # yes for anisotropic system
like membrane
useConstantRatio yes
langevin on
langevinDamping 5
langevinTemp $temp
langevinHydrogen off
# constant pressure
langevinPiston on
langevinPistonTarget 1.01325
langevinPistonPeriod 50.0
langevinPistonDecay 25.0
langevinPistonTemp $temp
run 100000
This archive was generated by hypermail 2.1.6 : Fri Dec 31 2021 - 23:17:08 CST