[gmx-users] Issue with PMF

Thu Apr 2 22:28:20 CEST 2020

On 4/2/20 4:19 PM, Alex wrote:
> Thanks Justin for the response.
>
> On Thu, Apr 2, 2020 at 3:21 PM Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>> On 4/2/20 12:53 PM, Alex wrote:
>>> Dear all,
>>> By using the below setting I am getting a very nice umbrella histogram
>> and
>>> iact also shows a very good integrated Autocorrelation time. However, the
>>> PMF profile is unexpected. Below I have also provided a gmx wham command
>> I
>>> am using, part of the pullx.xvg of one of the windows and also
>>> corresponding part in log file.
>>> Would you please help me find out what the problem might be?
>> What is the problem?
>>
> Below link shows the histogram and PMF, actually I was expecting that the
> PMF reach a plateau in the region marked by B in the figure as we are fully
> inside the bulk area of the thin film. However, there are some region from
> 15 to 20 kcal/mol different respect to the expected plateau. Here we have
> 10nm thick thin film.
>
> https://drive.google.com/open?id=1ve50thHG4wORe0-fKdoIsPncQpEyuwcz
>

Bootstrapping may underestimate the real error. How well converged is 
the PMF with respect to time, e.g. computing the PMF from the first half 
and last half of the simulation time?

>>> One more question: does it make sense to manually shift the plateau of
>> the
>>> PMF plot to zero (when we plot it), instead of using using the -zprof0
>> flag?
>>
>> There is no difference. The -zprof0 just adds a constant so the PMF is
>> zero at the specified point on the reaction coordinate.
>>
>>> %--------------------------------
>>> pull                     = yes
>>> pull-ngroups             = 2
>>> pull-ncoords             = 1
>>> pull-group1-name         = Mol_A
>>> pull-group2-name         = Thin_film
>>>
>>> pull-coord1-groups       = 1 2
>>> pull-coord1-type         = umbrella
>>> pull-coord1-dim          = N N Y
>>> pull-coord1-start        = no   ;--> manually we define the distance for
>>> all windows
>>> pull-coord1-rate         = 0.0  ;--> We don't want the roups to move
>>> pull-coord1-geometry     = distance
>>> pull-coord1-k            = 10000   ;;kJ/(mol nm^2)
>>> pull-print-components    = Yes
>>> pull-nstxout             = 2000
>>> pull-nstfout             = 2000
>>> pull-print-com1          = yes
>>> pull-coord1-init         = 1.650000
>>> %------------------------------------
>>> #Part of the one of the pullx.xvg file
>>> @ legend length 2
>>> @ s0 legend "1"
>>> @ s1 legend "1 dZ"
>>> @ s2 legend "1 g 1 Z"
>>> @ s3 legend "1 g 2 Z"
>>> 0.0000  0.0155171       0.0155171       9.89805 9.91357
>>> 2.0000  0.0103956       0.0103956       9.90227 9.91266
>>> 4.0000  0.0146962       0.0146962       9.89733 9.91202
>>> 6.0000  0.00656533      0.00656533      9.90404 9.91061
>>> 8.0000  0.00228074      0.00228074      9.90945 9.91174
>>> 10.0000 0.00967605      0.00967605      9.90514 9.91482
>>> ....
>>> %--------------------------------------------
>>> gmx_mpi wham -hist Histo.xvg -nBootstrap 200 -bins 200 -bs-method b-hist
>>> -bsres bsResult.xvg -bsprof bsProfs.xvg -if Fpull.dat -it TPR.dat -ac -o
>>> profile.xvg -b 3000 -unit kCal -v
>>>
>>>    %------------------------------------------------
>>> #Similar massage as below I have for all windows (tpr files).
>>> Reading file prd.1.tpr, VERSION 2018.7 (single precision)
>>>
>>> File prd.1.tpr, 1 coordinates, with these options:
>>>           Geometry distance   k = 10000     position = 0
>>   dimensions
>>> [N N Y] (1 dimensions). Used: Yes
>>>           Pull group coordinates of 2 groups expected in pullx files.
>>>           Reference value of the coordinate not expected in pullx files.
>>>
>>> Reading pull force file prd_pullf.1.xvg, expecting 2 columns:
>>>           Columns for pull coordinate 1
>>>                   reaction coordinate:             1
>>>                   center-of-mass of groups:        0
>>>                   reference position column:       No
>>>           Found 2501 times in prd_pullf.1.xvg
>>> %----------------------------------------------------------
> Actually, I was suspicious of these "Reference value of the coordinate not
> expected in pullx files."  and "reference position column:       No" which
> I have in log file of gmx wham for each windows's tpr.

I have no experience with that; the pullx.xvg file format has changed a 
lot over time. I always use pullf.xvg and it's never an issue.

> Concerning your PMF tutorial, when the distance between chain_A and chain_B
> is calculated by gmx distance, why the -oall is used to get the absolute
> distance out of x,y and z, while the -oxyz would be much better as we can
> take the Z component of the distance similar (corresponding) to the
> considered reaction coordinate in the tutorial?

The two quantities should be essentially equivalent for the purposes of 
the tutorial.

-Justin

> Regards,
> Alex
>
>>> Thank you,
>>> Alex
>> --
>> ==================================================
>>
>> Justin A. Lemkul, Ph.D.
>> Assistant Professor
>> Office: 301 Fralin Hall
>> Lab: 303 Engel Hall
>>
>> Virginia Tech Department of Biochemistry
>> 340 West Campus Dr.
>> Blacksburg, VA 24061
>>
>> jalemkul at vt.edu | (540) 231-3129
>> http://www.thelemkullab.com
>>
>> ==================================================
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-request at gromacs.org.
>>

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

==================================================