RE: pKa calculation with thermodynamic integration (TI)

From: Radak, Brian K (bradak_at_anl.gov)
Date: Wed May 24 2017 - 13:29:59 CDT

Sorry, those options are new in NAMD 2.12. They may be necessary depending on how you construct your thermodynamic cycle (or what you want to calculate).

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: Sadegh Faramarzi Ganjabad [safaramarziganjabad_at_mix.wvu.edu]
Sent: Wednesday, May 24, 2017 1:16 PM
To: Radak, Brian K
Cc: namd-l_at_ks.uiuc.edu; Hannes Loeffler
Subject: Re: namd-l: pKa calculation with thermodynamic integration (TI)

Brian,

I am using an example namd file from FEP tutorial. Now I put

alchElecLambdaStart 0.0

alchBondDecouple on
alchBondLambdaEnd 1.0

Interestingly, namd is complaining about the last two lines:

ERROR: The following variables were set in the
ERROR: configuration file but are NOT VALID
ERROR: alchBondLambdaEnd
ERROR: alchBondDecouple
FATAL ERROR: ERROR(S) IN THE CONFIGURATION FILE

I checked namd user guide, the syntax looks right. I have no idea why it gives error.

Sadegh

On Wed, May 24, 2017 at 12:35 PM, Radak, Brian K <bradak_at_anl.gov<mailto:bradak_at_anl.gov>> wrote:
Sadegh,

Without looking too carefully at your system, it's hard to say what is wrong. A few questions/comments:

1) did you mean to set "alchElecLambdaStart 1.0", should this not be 0.0?
2) you may also want to set "alchBondDecouple on" and "alchBondLambdaEnd 1.0" - this will turn on/off bonds between the atoms that are appearing/disappearing
3) does this work without pressure control?

Cheers,
Brian

Brian Radak
Postdoctoral Appointee
Leadership Computing Facility
Argonne National Laboratory

9700 South Cass Avenue, Bldg. 240
Argonne, IL 60439-4854
(630) 252-8643<tel:(630)%20252-8643>
brian.radak_at_anl.gov<mailto:brian.radak_at_anl.gov>
________________________________
From: Sadegh Faramarzi Ganjabad [safaramarziganjabad_at_mix.wvu.edu<mailto:safaramarziganjabad_at_mix.wvu.edu>]
Sent: Wednesday, May 24, 2017 10:48 AM
To: Radak, Brian K
Cc: namd-l_at_ks.uiuc.edu<mailto:namd-l_at_ks.uiuc.edu>; Hannes Loeffler

Subject: Re: namd-l: pKa calculation with thermodynamic integration (TI)

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<mailto: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<mailto: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<tel:%28630%29%20252-8643>
brian.radak_at_anl.gov<mailto:brian.radak_at_anl.gov>

________________________________________
From: owner-namd-l_at_ks.uiuc.edu<mailto:owner-namd-l_at_ks.uiuc.edu> [owner-namd-l_at_ks.uiuc.edu<mailto:owner-namd-l_at_ks.uiuc.edu>] on behalf of Hannes Loeffler [Hannes.Loeffler_at_stfc.ac.uk<mailto:Hannes.Loeffler_at_stfc.ac.uk>]
Sent: Wednesday, May 17, 2017 3:37 AM
To: namd-l_at_ks.uiuc.edu<mailto: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<tel:(630)%20252-8643>
brian.radak_at_anl.gov<mailto:brian.radak_at_anl.gov>

This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:21:18 CST