RE: Dealing with lone pairs

From: Pang, Yui Tik (andrewpang_at_gatech.edu)
Date: Wed Feb 03 2021 - 09:39:01 CST

Dear Seke,

Would you share your psf and pdb files, as well as the NAMD output log file so that we can better assess what the problem might be?

Best,
Andrew

From: Seke Keretsu<mailto:sekekeretsu_at_gmail.com>
Sent: Wednesday, February 3, 2021 4:05 AM
To: namd-l_at_ks.uiuc.edu<mailto:namd-l_at_ks.uiuc.edu>
Subject: namd-l: Dealing with lone pairs

Dear Experts,

I am posting this again as follow up to my previous post. Since then, I followed and tried all the advice but the error was not resolved.

I am using NAMD2.14 for a protein ligand complex in which there is a lone pair on the chlorobenzene of the ligand. The psf files were built in Charmm-Gui and solvated and ionized with VMD 1.9.4.

As suggested, I have added 'lonepairs on' to the configuration file.

In the psf file (ionized.psf), I added the lines for the one lone pair as advised (by Pang, Yui Tik):

         1 3 !NUMLP NUMLPH
         2 1 F 1.64000 0.00000 0.00000
      4737 4710 4697

Here in the second line, I have tried 1.64 and -1.64 for distance between lonepair and host. For polarized water it was negative, why?
In the third line, 4737 4710 4697 are the atomic ID of lonepair, host and atom connected to host.

I still get the error as shown below regardless of the temperature I use.

ERROR MESSAGE:

ENERGY: 50020 287.8329 875.2009 2562.5041 35.3830 -139424.9435 155009.5692 0.0000 0.0000 0.0000 19345.5467 0.0000 19345.5467 19345.5467 0.0000 123095.5120 107654.6219 396330.9410 123095.5120 107654.6219

TCL: Running for 500000 steps
ERROR: Margin is too small for 1 atom during timestep 50026.
ERROR: Incorrect nonbonded forces and energies may be calculated!
ERROR: Atom 4737 velocity is -2369.83 -18252.9 -15708.5 (limit is 12000, atom 395 of 507 on patch 37 pe 0)
ERROR: Atoms moving too fast; simulation has become unstable (1 atoms on patch 37 pe 0).
FATAL ERROR: Exiting prematurely; see error messages above.

As suggested (by Laurent Chaloin), I have also checked the content of "par_all36_cgenff.prm." Since I am using NAMD for simulation I assume the lines should not be comment-out.
LPH CLGR1 0.00 0.0000 ! aromatic halogen to lone pair
LPH BRGR1 0.00 0.0000 ! aromatic halogen to lone pair
LPH IGR1 0.00 0.0000 ! aromatic halogen to lone pair
However I could not find the following lines in my .prm file. I tried running a few times by appending the lines to my par_all36_cgenff.prm file. Did not change the outcome.

Please advise if I missed something here. Also, Is there an alternative to this, such as building a system without a lone pair without badly compromising my simulation protocol ?

Thank you.
Seke

CONFIGURATION FILE:

#############################################################
## JOB DESCRIPTION ##
#############################################################

# Minimization and Equilibration of
# Ubiquitin in a Water Box

#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################

structure ionized.psf
coordinates ionized.pdb

set temperature 100
set outputname equilibration02

firsttimestep 50000
bincoordinates equilibration.restart.coor
#binvelocities equilibration.restart.vel

#############################################################
## SIMULATION PARAMETERS ##
#############################################################

# Input
paraTypeCharmm on
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
temperature $temperature

# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
cutoff 12.0
switching on
switchdist 10.0
pairlistdist 14.0

# Integrator Parameters
timestep 1.0 ;# 2fs/step
rigidBonds all ;# needed for 2fs steps
nonbondedFreq 1
fullElectFrequency 2
stepspercycle 10

# Constant Temperature Control
langevin on ;# do langevin dynamics
langevinDamping 1 ;# damping coefficient (gamma) of 1/ps
langevinTemp $temperature
langevinHydrogen off ;# don't couple langevin bath to hydrogens

# Periodic Boundary Conditions

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

wrapAll on

# PME (for full-system periodic electrostatics)
PME yes
PMEGridSpacing 1.0

#manual grid definition
#PMEGridSizeX 45
#PMEGridSizeY 45
#PMEGridSizeZ 48

# Constant Pressure Control (variable volume)
useGroupPressure yes ;# needed for rigidBonds
useFlexibleCell no
useConstantArea no

langevinPiston on
langevinPistonTarget 1.01325 ;# in bar -> 1 atm
langevinPistonPeriod 100.0
langevinPistonDecay 50.0
langevinPistonTemp $temperature

constraints on
consexp 2
consref prot_posres.ref
conskfile prot_posres.ref
conskcol B
constraintScaling 10.0

# Output
outputName $outputname

restartfreq 500 ;# 500steps = every 1ps
dcdfreq 250
xstFreq 250
outputEnergies 100
outputPressure 100

lonepairs on

#############################################################
## EXTRA PARAMETERS ##
#############################################################

#############################################################
## EXECUTION SCRIPT ##
#############################################################

# Minimization
minimize 20 #system was minimized already
#reinitvels $temperature

run 500000 ;# 50ps

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