Re: pKa calculation with thermodynamic integration (TI)

From: Brian Radak (bradak_at_anl.gov)
Date: Tue May 30 2017 - 09:23:34 CDT

You probably due not need any pressure related flags is langevinPiston
is off.

As others asked previously (perhaps I derailed that line of inquiry),
can you verify that a single energy calculation ("run 0") gives a
sensible result? It is possible that your alchemical flags (the 1's and
-1's in the PDB) are not exactly set as you want them to be.

On 05/29/2017 03:45 PM, Sadegh Faramarzi Ganjabad wrote:
> Brian,
>
> I am using NAMD 2.12 now. I had two questions though
>
> 1. My system needs a minimization step before running free energy
> calculations. For FEP method there is a section for
> minimization+simulation in fep.tcl, namely runFEPmin. But there is not
> a similar section for TI in the same script.
>
> 2. In the configuration file I removed all lines in pressure control
> section and replaced them with the following
>
> useGroupPressure yes ;# needed for rigidBonds
> useFlexibleCell no
> useConstantArea no
>
> I thought that's all I need for an NVT system. Then I tried runFEPmin.
> Now I keep getting 'Segmentation fault' error that seems to be related
> to pressure of the system. Here are few lines from the output
>
> TCL: Original numsteps 1 will be ignored.
> TCL: Minimizing for 10000 steps
> 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
> TITITLE: TS BOND1 ELECT1 VDW1 BOND2
> ELECT2 VDW2
>
> ENERGY: 0 3109.7847 11048.2916 8201.0126
> 183.6201 -94866.8114 5152.1218 0.0000
> 0.0000 0.0000 -67171.9806 0.0000 -67171.9806
> -67171.9806 0.0000 -nan -nan
> 450474.3936 -nan -nan
> TI: 0 19.7158 7.3273 -8.2273 12.4099
> 22.8317 -8.3527
>
> OPENING EXTENDED SYSTEM TRAJECTORY FILE
> MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM
> LINE MINIMIZER REDUCING GRADIENT FROM nan TO nan
>
> Has anybody else had a similar issue?
>
> Thanks,
> Sadegh
>
> On Wed, May 24, 2017 at 2:29 PM, Radak, Brian K <bradak_at_anl.gov
> <mailto:bradak_at_anl.gov>> wrote:
>
> 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 <tel:%28630%29%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 1:16 PM
>
> *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,
>
> 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:%28630%29%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
> <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:%28630%29%20252-8643>
> brian.radak_at_anl.gov <mailto:brian.radak_at_anl.gov>
>
>
>
>

-- 
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:20 CST