Re: pKa calculation with thermodynamic integration (TI)

From: Sadegh Faramarzi Ganjabad (safaramarziganjabad_at_mix.wvu.edu)
Date: Wed May 24 2017 - 13:16:46 CDT

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> 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
> brian.radak_at_anl.gov
> ------------------------------
> *From:* Sadegh Faramarzi Ganjabad [safaramarziganjabad_at_mix.wvu.edu]
> *Sent:* Wednesday, May 24, 2017 10:48 AM
> *To:* Radak, Brian K
> *Cc:* 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> 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