Re: Doubts and problem in ABF simulations

From: Jerome Henin (heninj_at_gmail.com)
Date: Wed Sep 02 2009 - 15:34:54 CDT

Hi,
Just to rephrase what Giacomo said: a coordinate that works for SMD
might not work for ABF.
Jerome

On Wed, Sep 2, 2009 at 2:21 PM, Giacomo Fiorin<gfiorin_at_seas.upenn.edu> wrote:
> Hi Paulo, it definitely looks like, in some cases, we should allow the
> colvars module to handle the dummy atom with ABF as well.
>
> However, in your application: are you sure that the distance from the
> dummy atom is more appropriate?  Unless you fix the protein at an
> absolute position in space, the colvar you just defined will be the
> distance from the origin (or whatever point in space you set after
> dummyAtom), and NOT from the binding site of the protein.
>
> Also, take into account that when using SMD, you steer the ligand
> along a specific trajectory.  When using ABF over a distance variable,
> you also sample movements of the ligand orthogonal to that trajectory.
>  So the phase spaces covered with the two simulations will be
> different.
>
> Giacomo
>
>
> ---- ----
>  Giacomo Fiorin
>   ICMS - Institute for Computational Molecular Science
>     Temple University
>     1900 N 12 th Street, Philadelphia, PA 19122
>  work phone:   (+1)-215-204-4216
>  mobile:  (+1)-267-324-7676
>  mail:    giacomo.fiorin_at_gmail.com
> ---- ----
>
>
>
> On Wed, Sep 2, 2009 at 11:55 AM, Paulo Cesar Telles de
> Souza<paulocts_at_gmail.com> wrote:
>> Dear Jerome
>>
>> Thanks a lot for your help! I didn't know that the NAMD 2.7b1 altered ABF
>> method (and others) in such a way. I have decided to use this new version of
>> NAMD with colvars module. I don't think that it is necessary employing
>> multidimensional ABF for my problem. I have obtained PMF using a SMDs and
>> Jarzynski equality and now I would like to compare the result with one
>> obtained using another free energy calculation method. To specify the same
>> reaction coordinate that I have used in SMD, I have tried to use the
>> distance component, that is, I chose ligand atoms for group 1 and dummyAtom
>> option for group 2 (setting the appropriate coordinate to reproduce the same
>> reaction coordinate), but, still, it didn't work. I pasted the error message
>> below. Is it necessary to use another option  together with dummyAtom?
>>
>> Bests,
>> Paulo
>>
>> colvars:
>> ----------------------------------------------------------------------
>> colvars: Initializing the collective variables module, version 29-09-2008.
>> colvars: # colvarsTrajFrequency = 500
>> colvars: # colvarsRestartFrequency = 500
>> colvars: # trajAppend = off [default]
>> colvars: The restart output state file will be "protabf-a.colvars.state".
>> colvars: The trajectory file will be "protabf-a.colvars.traj".
>> colvars: The final output state file will be "protabf-a.colvars.state".
>> colvars: # analysis = off [default]
>> colvars:
>> ----------------------------------------------------------------------
>> colvars:   Initializing a new collective variable.
>> colvars:   # name = Distance
>> colvars:   Initializing a new "distance" component.
>> colvars:     # componentCoeff = 1 [default]
>> colvars:     # componentExp = 1 [default]
>> colvars:     # oneSiteSystemForce = off [default]
>> colvars:       Initializing atom group "group1".
>> colvars:       Atom group "group1" defined, 35 atoms, total mass = 650.977.
>> colvars:       Initializing atom group "group2".
>> colvars:       Atom group "group2" defined, 0 atoms, total mass = 1.
>> colvars:   All components initialized.
>> colvars:   # width = 0.1
>> colvars:   # lowerBoundary = 16
>> colvars:   Lower boundary defined.
>> colvars:   # upperBoundary = 31
>> colvars:   Upper boundary defined.
>> colvars:   # lowerWallConstant = 10
>> colvars:   Applying a harmonic lower wall at 16, with coefficient 10.
>> colvars:   # upperWallConstant = 10
>> colvars:   Applying a harmonic upper wall at 31, with coefficient 10.
>> colvars:   # extendedLagrangian = off [default]
>> colvars:   # outputValue = on [default]
>> colvars:   # outputVelocity = off [default]
>> colvars:   # outputSystemForce = off [default]
>> colvars:   # outputAppliedForce = off [default]
>> colvars: ------------------------------
>> ----------------------------------------
>> colvars: Collective variables initialized, 1 in total.
>> colvars:
>> ----------------------------------------------------------------------
>> colvars:   Initializing a new "abf" instance.
>> colvars:   # name = "abf1" [default]
>> colvars:   # colvars = { Distance }
>> colvars:   # applybias = on [default]
>> colvars:   # hidejacobian = off [default]
>> colvars:   Jacobian (geometric) forces will be included in reported free
>> energy gradients.
>> colvars:   # fullsamples = 500
>> colvars:   # inputprefix =  [default]
>> colvars:   # outputfreq = 500 [default]
>> colvars:   Finished ABF setup.
>> colvars:
>> ----------------------------------------------------------------------
>> colvars: Collective variables biases initialized, 1 in total.
>> colvars:
>> ----------------------------------------------------------------------
>> colvars: Collective variables module initialized.
>> colvars:
>> ----------------------------------------------------------------------
>> Info: Startup phase 7 took 0.012846 s, 7.72452 MB of memory in use
>> Info: Startup phase 8 took 0.000380039 s, 12.3929 MB of memory in use
>> Info: Finished startup at 0.541487 s, 12.3929 MB of memory in use
>>
>> TCL: Running for 4000000 steps
>> ETITLE:      TS           BOND          ANGLE          DIHED
>> IMPRP               ELECT            VDW       BOUNDARY
>> MISC        KINETIC               TOTAL           TEMP
>> POTENTIAL         TOTAL3        TEMPAVG            PRESSURE
>> GPRESSURE         VOLUME       PRESSAVG      GPRESSAVG
>>
>> ENERGY:       0       849.8375      2408.4881      1348.5871
>> 135.3159          -5237.1600     -1091.2560         0.0000
>> 0.0000      3104.3467           1518.1595       292.6250     -1586.1873
>> 1557.9272       292.6250            112.1152       111.0410
>> 590148.0000       112.1152       111.0410
>>
>> colvars:   Error: system forces are not available from a dummy atom group.
>> colvars:   If this error message is unclear, try recompile with
>> -DCOLVARS_DEBUG.
>> FATAL ERROR: Error in the collective variables module: exiting.
>>
>> ------------- Processor 0 Exiting: Called CmiAbort ------------
>> Reason: FATAL ERROR: Error in the collective variables module: exiting.
>>
>>
>> [0] Stack Traceback:
>>   [0] CmiAbort+0x7f  [0xa4526d]
>>   [1] _Z8NAMD_diePKc+0x62  [0x500082]
>>   [2] _ZN16colvarproxy_namd11fatal_errorERKSs+0x64  [0x844b68]
>>   [3] _ZN12colvarmodule11fatal_errorERKSs+0x18  [0x79ae22]
>>   [4] _ZNK12colvarmodule10atom_group12system_forceEv+0x53  [0x816f75]
>>   [5] _ZN6colvar8distance19calc_force_invgradsEv+0xd7  [0x829145]
>>   [6] _ZN6colvar4calcEv+0xc03  [0x7cca1f]
>>   [7] _ZN12colvarmodule4calcEv+0x6c  [0x7969bc]
>>   [8] _ZN16colvarproxy_namd9calculateEv+0x497  [0x8415a1]
>>   [9] _ZThn16_N16colvarproxy_namd9calculateEv+0xa  [0x841108]
>>   [10]
>> _ZN12GlobalMaster11processDataEPiS0_P6VectorS2_S2_S0_S0_S2_S0_S0_S2_+0x71
>> [0x77b165]
>>   [11] _ZN18GlobalMasterServer11callClientsEv+0x469  [0x77f191]
>>   [12] _ZN18GlobalMasterServer8recvDataEP20ComputeGlobalDataMsg+0x6f7
>> [0x77e7d7]
>>   [13] _ZN10ComputeMgr21recvComputeGlobalDataEP20ComputeGlobalDataMsg+0x12
>> [0x569388]
>>   [14]
>> _ZN18CkIndex_ComputeMgr48_call_recvComputeGlobalData_ComputeGlobalDataMsgEPvP10ComputeMgr+0xf
>> [0x569373]
>>   [15] CkDeliverMessageFree+0x21  [0x9c2d71]
>>   [16] _Z15_processHandlerPvP11CkCoreState+0x509  [0x9c2365]
>>   [17] CsdScheduleForever+0xa5  [0xa4bcf5]
>>   [18] CsdScheduler+0x1c  [0xa4b8f6]
>>   [19] _ZN7BackEnd7suspendEv+0xb  [0x508add]
>>   [20] _ZN9ScriptTcl7Tcl_runEPvP10Tcl_InterpiPPc+0x140  [0x8fcde4]
>>   [21] TclInvokeStringCommand+0x91  [0xa6e208]
>>   [22] namd2 [0xaa4058]
>>   [23] Tcl_EvalEx+0x176  [0xaa469b]
>>   [24] Tcl_EvalFile+0x134  [0xa9c0a4]
>>   [25] _ZN9ScriptTcl3runEPc+0x14  [0x8fc4e2]
>>   [26] _Z18after_backend_initiPPc+0x22b  [0x50481b]
>>   [27] main+0x3a  [0x5045ba]
>>   [28] __libc_start_main+0xe6  [0x7ff4b89e8466]
>>   [29] _ZNSt8ios_base4InitD1Ev+0x52  [0x4ff9ea]
>> Charm++ fatal error:
>> FATAL ERROR: Error in the collective variables module: exiting.
>>
>>
>> Aborted
>> 2009/8/31 Jerome Henin <jhenin_at_cmm.chem.upenn.edu>
>>>
>>> Dear Paulo,
>>>
>>> About your doubt: it is not easy to describe ligand-receptor binding
>>> with a single parameter. Sop far, I am not aware of anyone
>>> successfully applying ABF to a protein-ligand binding problem. You may
>>> have better luck with a multidimensional calculation, using either ABF
>>> or metadynamics.
>>>
>>> About your problem: indeed, there seems to be an error in the
>>> documentation for the abscissa parameter. As mentioned prominently on
>>> the webpage you link to, the NAMD 2.6 ABF code is deprecated, and only
>>> the colvars module of NAMD 2.7b1 is actively maintained. If you cannot
>>> use NAMD 2.7b1, one option is to modify the abscissa.tcl script to
>>> accept an atom group for abf3 rather than a single atom. This should
>>> be relatively straightforward. However, it will not solve the problem
>>> of choosing the right coordinate(s) for the job.
>>>
>>> Best,
>>> Jerome
>>>
>>>
>>> On Mon, Aug 31, 2009 at 12:00 PM, Paulo Cesar Telles de
>>> Souza<paulocts_at_gmail.com> wrote:
>>> > Dear all,
>>> >
>>> > I have some doubts and one problem on running ABF simulations in NAMD.
>>> >
>>> > Doubts: I'm not sure what order parameter is better for my system,
>>> > mainly
>>> > because I have tested to use the z-coord order parameter and it didn't
>>> > work.
>>> > I selected some atoms in the binding site (abf1) and the ligant atoms
>>> > (abf2). The z axis of the box is parallel to the direction from the
>>> > ligand
>>> > in the receptor to the solvent bulk and the distance vector between abf1
>>> > and
>>> > abf2 is almost parallel to the z-direction. For my suprise, when I
>>> > performed
>>> > the ABF simulation, the ligant "walk'' in the protein surface. Did I
>>> > choose
>>> > the wrong order parameter or the problem is due to my vector abf1-abf2 ?
>>> > Does the vector abf1-abf2 change during ABF simulations ?
>>> >
>>> > Problem: I have tried to use the abscissa order parameter to obtain a
>>> > PMF
>>> > along a dissociation pathway through which a ligand unbinds from its
>>> > receptor (the reaction coordinate is ligant-protein to ligant-solvent).
>>> > I
>>> > used a group of atoms in abf3 but received the following message error.
>>> > This
>>> > is strange because the manual and this repository
>>> > (http://www.edam.uhp-nancy.fr/ABF/repository.html) specify that I can
>>> > choose
>>> > a list of atoms. If I try only with one atom, it works. How can I use a
>>> > list
>>> > in abf3 ( I want that a entire system to move along the vector defined
>>> > by
>>> > abf1 and abf2) ?
>>> >
>>> > Thanks on advance
>>> >
>>> > Paulo
>>> >
>>> > TCL: ABF> ---------------------------------------------
>>> > TCL: ABF> Adaptive Biasing Force protocol version 1.7
>>> > TCL: ABF> ---------------------------------------------
>>> > TCL: ABF>
>>> > TCL: ABF> Using coordinate type : abscissa
>>> > TCL: ABF> Position of an atom along the axis joining two groups
>>> > TCL: ABF>              dxi : 0.1
>>> > TCL: ABF>          outFile : abf.dat        [default]
>>> > TCL: ABF>      fullSamples : 500
>>> > TCL: ABF>     writeFxiFreq : 0              [default]
>>> > TCL: ABF>      writeXiFreq : 100
>>> > TCL: ABF>             abf2 : 4240 4241 4242 4243 4244 4245
>>> > TCL: ABF>             fMax : 60.0           [default]
>>> > TCL: ABF>       outputFreq : 5000           [default]
>>> > TCL: ABF>          inFiles :                [default]
>>> > TCL: ABF>      historyFile : none           [default]
>>> > TCL: ABF>               df : 1.0            [default]
>>> > TCL: ABF>     moveBoundary : 0              [default]
>>> > TCL: ABF>         distFile : none           [default]
>>> > TCL: ABF>             temp : 300.0          [default]
>>> > TCL: ABF>       forceConst : 10.0
>>> > TCL: ABF>          dSmooth : 0.1
>>> > TCL: ABF>            xiMin : -1.0
>>> > TCL: ABF>            xiMax : 15.0
>>> > TCL: ABF>             abf1 : 2724 2725 2726 2727 2728 2729 2730 2731
>>> > 2732
>>> > 2733 2734 2735 2736 2737 2738 2739 2806 2807 2808 2809 2810 2811 2812
>>> > 2813
>>> > 2814 2815 2816 2817 2818 2819 2820 2821 2822 2823 2824 2825 2826 2827
>>> > 2828
>>> > 2829 2830 2831 3040 3041 3042 3043 3044 3045 3046 3047 3048 3049 3050
>>> > 3051
>>> > 3052 3053 3054 3055 3056 3103 3104 3105 3106 3107 3108 3109 3110 3111
>>> > 3112
>>> > 3113 3114 3115 3116 3117 3118 3119 3120 3121 3418 3419 3420 3421 3422
>>> > 3423
>>> > 3424 3425 3426 3427 3428 3429 3430 3431 3432 3433 3434 3435 3436 3437
>>> > 3438
>>> > 3439 3440 3441 3497 3498 3499 3500 3501 3502 3503 3504 3505 3506 3507
>>> > 3508
>>> > 3509 3510 3511 3512 3513 3514 3515
>>> > TCL: ABF> Accumulating force data in 160 bins
>>> > TCL: ABF> Moving abf.dat to abf.dat.BAK
>>> > TCL: ABF> Using the abscissa of atom 4240 4243
>>> > TCL: expected integer but got "4240 4243"
>>> > FATAL ERROR: expected integer but got "4240 4243"
>>> >     while executing
>>> > "addatom $abf3"
>>> >     (in namespace eval "::ABF::ABFcoord" script line 9)
>>> >     invoked from within
>>> > "namespace eval ABFcoord {
>>> >
>>> >     set abf1 $::ABF::abf1
>>> >     set abf2 $::ABF::abf2
>>> >     set abf3 $::ABF::abf3
>>> >
>>> >     print "ABF> Using the abscissa of atom $abf3"
>>> >
>>> >     adda..."
>>> >     (procedure "ABFstartup" line 2)
>>> >     invoked from within
>>> > "ABFstartup"
>>> >     (in namespace eval "::ABF" script line 4)
>>> >     invoked from within
>>> > "namespace eval ::ABF {
>>> >
>>> > # coordinate-specific startup procedure
>>> > ABFstartup
>>> >
>>> > if { [info exists restraintList] } {
>>> >
>>> >     array set rArray $restraintList
>>> > ..."
>>> >     (file "/home/paulo/Desktop/Tutorial-ABF/abf-1.8/abf_script.tcl" line
>>> > 766)
>>> > ------------- Processor 0 Exiting: Called CmiAbort ------------
>>> > Reason: FATAL ERROR: expected integer but got "4240 4243"
>>> >     while executing
>>> > "addatom $abf3"
>>> >     (in namespace eval "::ABF::ABFcoord" script line 9)
>>> >     invoked from within
>>> > "namespace eval ABFcoord {
>>> >
>>> >     set abf1 $::ABF::abf1
>>> >     set abf2 $::ABF::abf2
>>> >     set abf3 $::ABF::abf3
>>> >
>>> >     print "ABF> Using the abscissa of atom $abf3"
>>> >
>>> >     adda..."
>>> >     (procedure "ABFstartup" line 2)
>>> >     invoked from within
>>> > "ABFstartup"
>>> >     (in namespace eval "::ABF" script line 4)
>>> >     invoked from within
>>> > "namespace eval ::ABF {
>>> >
>>> > # coordinate-specific startup procedure
>>> > ABFstartup
>>> >
>>> > if { [info exists restraintList] } {
>>> >
>>> >     array set rArray $restraintList
>>> > ..."
>>> >     (file "/home/paulo/Desktop/Tutorial-ABF/abf-1.8/abf_script.tcl" line
>>> > 766)
>>> >
>>> > [0] Stack Traceback:
>>> >   [0] CmiAbort+0x7f  [0xa4526d]
>>> >   [1] _Z8NAMD_diePKc+0x62  [0x500082]
>>> >   [2] _ZN15GlobalMasterTcl10initializeEv+0x130  [0x789f22]
>>> >   [3] _ZN15GlobalMasterTclC9Ev+0x66  [0x786426]
>>> >   [4] _ZN15GlobalMasterTclC1Ev+0x6  [0x78646c]
>>> >   [5] _ZN10ComputeMgr14createComputesEP10ComputeMap+0x103  [0x56cb59]
>>> >   [6] _ZN4Node7startupEv+0x287  [0x89f303]
>>> >   [7] _ZN12CkIndex_Node18_call_startup_voidEPvP4Node+0x12  [0x89f078]
>>> >   [8] CkDeliverMessageFree+0x21  [0x9c2d71]
>>> >   [9] _Z15_processHandlerPvP11CkCoreState+0x509  [0x9c2365]
>>> >   [10] CsdScheduleForever+0xa5  [0xa4bcf5]
>>> >   [11] CsdScheduler+0x1c  [0xa4b8f6]
>>> >   [12] _ZN7BackEnd7suspendEv+0xb  [0x508add]
>>> >   [13] _ZN9ScriptTcl9initcheckEv+0x80  [0x90041e]
>>> >   [14] _ZN9ScriptTcl7Tcl_runEPvP10Tcl_InterpiPPc+0x20  [0x8fccc4]
>>> >   [15] TclInvokeStringCommand+0x91  [0xa6e208]
>>> >   [16] /home/paulo/programas/NAMD/namd2 [0xaa4058]
>>> >   [17] Tcl_EvalEx+0x176  [0xaa469b]
>>> >   [18] Tcl_EvalFile+0x134  [0xa9c0a4]
>>> >   [19] _ZN9ScriptTcl3runEPc+0x14  [0x8fc4e2]
>>> >   [20] _Z18after_backend_initiPPc+0x22b  [0x50481b]
>>> >   [21] main+0x3a  [0x5045ba]
>>> >   [22] __libc_start_main+0xe6  [0x7ff1a10d4466]
>>> >   [23] _ZNSt8ios_base4InitD1Ev+0x52  [0x4ff9ea]
>>> > Fatal error on PE 0> FATAL ERROR: expected integer but got "4240 4243"
>>> >     while executing
>>> > "addatom $abf3"
>>> >     (in namespace eval "::ABF::ABFcoord" script line 9)
>>> >     invoked from within
>>> > "namespace eval ABFcoord {
>>> >
>>> >     set abf1 $::ABF::abf1
>>> >     set abf2 $::ABF::abf2
>>> >     set abf3 $::ABF::abf3
>>> >
>>> >     print "ABF> Using the abscissa of atom $abf3"
>>> >
>>> >     adda..."
>>> >     (procedure "ABFstartup" line 2)
>>> >     invoked from within
>>> > "ABFstartup"
>>> >     (in namespace eval "::ABF" script line 4)
>>> >     invoked from within
>>> > "namespace eval ::ABF {
>>> >
>>> > # coordinate-specific startup procedure
>>> > ABFstartup
>>> >
>>> > if { [info exists restraintList] } {
>>> >
>>> >     array set rArray $restraintList
>>> > ..."
>>> >     (file "/home/paulo/Desktop/Tutorial-ABF/abf-1.8/abf_script.tcl" line
>>> > 766)
>>
>>
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:53:15 CST