Re: Dealing with lone pairs

From: David Hardy (dhardy_at_ks.uiuc.edu)
Date: Wed Feb 03 2021 - 14:31:29 CST

Hi Seke,

Please at least post the log file (redirected output from NAMD run).

Even better if you could post all of the simulation files so that we can reproduce your error.

The error message seems to indicate that NAMD is not actually treating the particle as a lone pair, but we can’t tell without seeing the actual output.

Thanks,
Dave

> On Feb 3, 2021, at 2:51 AM, Seke Keretsu <sekekeretsu_at_gmail.com> wrote:
>
> 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