Stability issues with POPC membrane using Drude force field

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