Re: pKa calculation with thermodynamic integration (TI)

From: Sadegh Faramarzi Ganjabad (safaramarziganjabad_at_mix.wvu.edu)
Date: Tue May 30 2017 - 13:12:23 CDT

Brian,

I get the same error when I run it for 0 timesteps. Then I tried minimizing
it first and then running TI. It gives the following error

FATAL ERROR: Low global exclusion count! (125865 vs 125915) System
unstable or pairlistdist or cutoff too small.

I hope it's not because of the dual topology. I am going to post it in a
new thread. My system might be out of equilibrium, regardless of the TI set
up.

Thanks,
Sadegh

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

> 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> 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 <%28630%29%20252-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 <%28630%29%20252-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 <%28630%29%20252-8643>
>>>> 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 : Sun Dec 31 2017 - 23:21:19 CST