Re: pKa calculation with thermodynamic integration (TI)

From: Sadegh Faramarzi Ganjabad (safaramarziganjabad_at_mix.wvu.edu)
Date: Wed May 24 2017 - 10:48:29 CDT

Brian,

Sorry I forgot the parameter file you sent me.

When I try to run TI, I get the following error

FATAL ERROR: Low global exclusion count! (125820 vs 125908) System
unstable or pairlistdist or cutoff too small. Presumably it's related to
pressure of the system. Here is a part of my log file

TI: WINDOW TO HAVE 100 STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.
PRESSURE: 0 -nan -nan -nan -nan -nan -nan -nan -nan -nan
GPRESSURE: 0 -nan -nan -nan -nan -nan -nan -nan -nan -nan
ETITLE: TS BOND ANGLE DIHED IMPRP
              ELECT VDW BOUNDARY MISC
 KINETIC TOTAL TEMP POTENTIAL TOTAL3
     TEMPAVG PRESSURE GPRESSURE VOLUME
PRESSAVG GPRESSAVG

ENERGY: 0 3117.7968 11053.5678 8207.2710 183.7892
        -94787.6786 5151.4818 0.0000 0.0000
-nan -nan -nan -67073.7721 nan
    -nan -nan -nan 450474.3936 -nan
        -nan

OPENING TI ENERGY OUTPUT FILE
OPENING EXTENDED SYSTEM TRAJECTORY FILE
FATAL ERROR: Low global exclusion count! (125820 vs 125908) System
unstable or pairlistdist or cutoff too small.

FATAL ERROR: See http://www.ks.uiuc.edu/Research/namd/bugreport.html

And here is my namd configuration file

# MD SETUP

timestep 2.0
numsteps 1

# FLEXIBLE CELL

useflexiblecell no

set temperature 310

# INPUT

structure all-ion-beta.psf
parameters par_all36m_prot.prm
parameters par_cph36_prot.prm
parameters par_all36_cgenff.prm
parameters par_all36_na.prm
parameters par_all36_carb.prm
parameters par_all36_lipid.prm
parameters toppar_water_ions_fixed.str
parameters retinal-pro.str
#parameters toppar_all36_lipid_model.str

paraTypeCharmm on
temperature $temperature

coordinates all-ion-beta.pdb
#bincoordinates equilibration.coor
#binvelocities equilibration.vel

#extendedSystem equilibration.xsc

# OUTPUT

outputenergies 100
outputtiming 100
outputpressure 100
binaryoutput yes

outputname forward-noshift-pr
restartname forward-noshift-pr
restartfreq 100
binaryrestart yes

XSTFreq 100

# PME

PME yes
PMETolerance 10e-6
PMEInterpOrder 4

PMEGridSpacing 1.0

# Periodic Boundary Conditions
if {1} {
cellBasisVector1 73.60 0.0 0.0
cellBasisVector2 0.0 73.60 0.0
cellBasisVector3 0.0 0 83.16
cellOrigin 0.030350102111697197 -0.059581365436315536
0.2382895052433014
}

# WRAP
wrapAll on

# CONSTANT-T

langevin on
langevinTemp 310.0
langevinDamping 1.0

# CONSTANT-P

LangevinPiston on
LangevinPistonTarget 1
LangevinPistonPeriod 75
LangevinPistonDecay 25
LangevinPistonTemp 300
StrainRate 0.0 0.0 0.0
useGroupPressure yes

# SPACE PARTITIONING

splitpatch hydrogen
hgroupcutoff 2.8
stepspercycle 20
margin 1.0

# CUT-OFFS

switching on
switchdist 10.0
cutoff 12.0
pairlistdist 12.0

# RESPA

fullElectFrequency 2
nonbondedFreq 1

# 1-4 NON-BONDED

exclude scaled1-4
1-4scaling 1.0

# COM

commotion no

# SHAKE

rigidbonds all
rigidtolerance 0.000001
rigiditerations 400

# EXTRA BONDS

extraBonds on
extraBondsFile ref
########################

# FEP PARAMETERS

source ../tools/fep.tcl

alch on
alchType TI
alchFile all-ion-beta.pdb
alchCol B
alchOutFile forward-noshift-pr.fepout
alchOutFreq 10

alchVdwLambdaEnd 1.0
alchElecLambdaStart 1.0
alchVdWShiftCoeff 0.0
alchDecouple no

alchEquilSteps 100
set numSteps 500

runFEP 0.0 1.0 0.0625 $numSteps

#######################

It happens even when I switch LangevinPiston off. I am applying harmonic
restraints with extraBonds. I am wondering how I can deal with this error.
Any help is appreciated.

Thanks,
Sadegh

On Tue, May 23, 2017 at 10:32 AM, Brian Radak <bradak_at_anl.gov> wrote:

> Those are defined in the parameter file that I sent (par_cph36_prot.prm).
> In principle, many of those parameters do not affect the final value, but
> DO affect the rate of convergence. Just to be clear, I cannot, at this
> time, offer any guarantees on the latter - they are experimental.
> The only important thing is that the dummy atoms:
> 1) do not have Lennard-Jones interactions with other particles
> 2) only have minimal bonded interactions with the rest of the system
>
> Both of those requirements are met by the parameters I defined.
>
> Cheers,
> Brian
>
>
>
> On 05/23/2017 01:04 AM, Sadegh Faramarzi Ganjabad wrote:
>
> Brian,
>
> I already used "set beta" keyword of VMD to assign values for beta. Now
> I'm wondering about parameters for dummy atoms. NAMD complains about not
> finding vdw parameter for HD atom type, for example. Do I have to manually
> assign a value (say 0) for that?
>
> Thanks,
> Sadegh
>
> On Wed, May 17, 2017 at 10:54 AM, Radak, Brian K <bradak_at_anl.gov> wrote:
>
>> You can do all of this from inside psfgen now:
>>
>> patch GLUHD <segid:resid>
>> regenerate angles dihedrals
>> guesscoord
>>
>> set atom0List {CG HG1 HG2 CD OE1 HE1 OE2 HE2}
>> foreach atom0 $atom0List {
>> set atom1 [format "%d1" $atom0] ; # this only works bc of the arbitrary
>> naming convention in the patch
>> psfset beta $atom0 -1.0
>> psfset beta $atom1 1.0
>> }
>>
>> writepsf <psfname>
>> writepdb <pdbname> ;# this assigns atom indices to all of the atoms
>>
>> set ebondsfile [open "foo.extrabonds" "w"]
>> foreach atom0 $atom0List {
>> set atom1 [format "%d1" $atom0] ;
>> set i [expr {segment atomid <segid> <resid> $atom0] - 1}] ;# note the
>> index shift
>> set j [expr {segment atomid <segid> <resid> $atom1] - 1}]
>> puts $ebondsfile [format "bond %d %d %f %f" $i $j 100.0 0.0]
>> }
>> close $ebondsfile
>>
>> Brian Radak
>> Postdoctoral Appointee
>> Leadership Computing Facility
>> Argonne National Laboratory
>>
>> 9700 South Cass Avenue, Bldg. 240
>> Argonne, IL 60439-4854
>> (630) 252-8643
>> brian.radak_at_anl.gov
>>
>> ________________________________________
>> From: owner-namd-l_at_ks.uiuc.edu [owner-namd-l_at_ks.uiuc.edu] on behalf of
>> Hannes Loeffler [Hannes.Loeffler_at_stfc.ac.uk]
>> Sent: Wednesday, May 17, 2017 3:37 AM
>> To: namd-l_at_ks.uiuc.edu; Sadegh Faramarzi Ganjabad
>> Subject: Re: namd-l: pKa calculation with thermodynamic integration (TI)
>>
>> > Is that what I am supposed to do? if so, how can I tight (for
>> > example) CD and CD1 atoms together by a harmonic restraint? at the
>> > beginning, there is no distance between them.
>>
>> You need to mark the outgoing group with -1 in the column that you have
>> chosen (alchCol). Look into the extraBonds feature to set additional
>> restraints.
>>
>>
>> Cheers,
>> Hannes.
>>
>>
>
> --
> Brian Radak
> Postdoctoral Appointee
> Leadership Computing Facility
> Argonne National Laboratory
>
> 9700 South Cass Avenue, Bldg. 240
> Argonne, IL 60439-4854
> (630) 252-8643
> brian.radak_at_anl.gov
>

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:20:19 CST