Dealing with lone pairs (LP) during simulation

From: Seke Keretsu (
Date: Sat Jan 30 2021 - 01:05:09 CST

Dear Expert,

I came across the following error during the equilibration of my protein
ligand system which has a CL atom (in ligand). I prepared the system using

ERROR: Atom 4737 velocity is -11982.1 12908.1 -3406.91 (limit is 12000,
atom 983 of 1008 on patch 17 pe 1)
ERROR: Atoms moving too fast; simulation has become unstable (1 atoms on
patch 17 pe 1).

Here, the atom no. 4737 corresponds to the lone pair connected to CL atom
in the system.

I have checked the mailing list which suggested to perform the simulation
gradually from low temperature for atoms moving fast. However, I could not
find a solution for dealing with lone pairs. I keep getting the same error
regardless of the temperature.

Is there a command I need to give in the configuration file?
Please give me some advice or link to a solution.

The configuration file used during equilibration is given below.
structure ionized.psf
coordinates ionized.pdb

set temp 50;
set outputname equilibration03;

#source step5_charmm2namd.str
set inputname equilibration;
outputname $outputname;
binCoordinates equilibration.restart.coor;
binVelocities equilibration.restart.vel;
extendedSystem equilibration.restart.xsc;

firsttimestep 50000;
restartfreq 1000; # 1000 steps = every 2ps
dcdfreq 5000;
dcdUnitCell yes; # the file will contain unit
cell info in the style of
                                            # charmm dcd files. if yes, the
dcd files will contain
                                            # unit cell information in the
style of charmm DCD files.
xstFreq 5000; # XSTFreq: control how often
the extended systen configuration
                                            # will be appended to the XST
outputEnergies 125; # 125 steps = every 0.25ps
                                            # The number of timesteps
between each energy output of NAMD
outputTiming 1000; # The number of timesteps
between each timing output shows
                                            # time per step and time to

# Force-Field Parameters
paraTypeCharmm on; # We're using charmm type
parameter file(s)
                                            # multiple definitions may be
used but only one file per definition
parameters jz4/jz4.prm
parameters toppar/par_all36_carb.prm
parameters toppar/par_all36_cgenff.prm
parameters toppar/par_all36_lipid.prm
parameters toppar/par_all36m_prot.prm
parameters toppar/par_all36_na.prm
parameters toppar/par_interface.prm
parameters toppar_water_ions.str
# 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 yes; # New option for force-based
switching of vdW
                                            # if both switching and
vdwForceSwitching are on CHARMM force
                                            # switching is used for vdW

# 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 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
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
rigidBonds all; # Bound constraint all bonds
involving H are fixed in length
nonbondedFreq 1; # nonbonded forces every step
fullElectFrequency 1; # PME every step

# Constant Temperature Control ONLY DURING EQUILB
reassignFreq 500; # reassignFreq: use this to
reassign velocity every 500 steps
reassignTemp $temp;

cellBasisVector1 69.61299896240234 0 0
cellBasisVector2 0 76.30400085449219 0
cellBasisVector3 0 0 74.61400032043457
cellOrigin 12.257286071777344 7.4436726570129395 -19.06593894958496

wrapWater on; # wrap water to central cell
wrapAll on; # wrap other molecules too
#wrapNearest $wrapnearst; # 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; # keeps the ratio of the unit
cell in the x-y plane constant A=B

langevin on
langevinDamping 1.0
langevinTemp $temp
langevinHydrogen off

constraints on
consexp 2
consref prot_posres.ref
conskfile prot_posres.ref
conskcol B
constraintScaling 5.0
minimize 5000
run 125000

Thank you.

Best, Seke

This archive was generated by hypermail 2.1.6 : Fri Dec 31 2021 - 23:17:10 CST