Re: pKa calculation with thermodynamic integration (TI)

From: Sadegh Faramarzi Ganjabad (safaramarziganjabad_at_mix.wvu.edu)
Date: Mon May 29 2017 - 15:45:10 CDT

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> 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
> 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> 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 : Sun Dec 31 2017 - 23:21:19 CST