[gmx-users] gmx wham (again)

Alex nedomacho at gmail.com
Sat Apr 21 03:39:01 CEST 2018 

It is definitely not (confinement of the pore on one hand, bulk solvent on
the other), but WHAM bases all of its calculations on the forces or
displacements in the orthogonal direction... I do suspect there's something
wrong with using the restraints here, but aside from some form of coupling
between components can't put my finger on it... How this coupling would
arise, and in what amount, I don't really have a good feeling for that. And
that magic cylinder isn't working. What to do?.. :)

On Fri, Apr 20, 2018 at 7:25 PM, Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> On 4/19/18 11:25 PM, Alex wrote:
>
>> Justin,
>>
>> Since that cylinder thing is buggy and I still haven't been able to get a
>> stable run (and I recall this is why I decided not to use it earlier), back
>> to the restraints...
>> If the WHAM procedure produces an essential relative energy fit and the
>> radial restraints are: 1) present in every config, including bulk, 2) in
>> the direction orthogonal to the reaction coordinate, and 3) 10 times lower
>> than e.g. the force constant corresponding to the ion-pore binding I am
>> after, is there much to account for in terms of artifacts? This is
>> basically a k<r^2>/2 energy term present in every config, where r is
>> perpendicular to the reaction coordinate. Anything I am missing in terms of
>> generality?
>>
>>
> You're assuming that the restraint does the same amount of work in every
> window, which may not be the case.
>
> -Justin
>
>
> Thank you,
>>
>> Alex
>>
>> On 4/19/2018 5:28 AM, Justin Lemkul wrote:
>>
>>>
>>>
>>> On 4/18/18 6:30 PM, Alex wrote:
>>>
>>>> An otherwise perfectly well-behaved setup is doing something
>>>> catastrophic.
>>>> The pull code is:
>>>>
>>>> ; Pull code
>>>> pull                 = yes
>>>> pull-coord1-type     = umbrella
>>>> pull_ngroups            = 2
>>>> pull_ncoords            = 1
>>>> pull_group1_name        = K
>>>> pull_group2_name        = CNT
>>>> pull_coord1_geometry    = cylinder
>>>> pull-cylinder-r         = 1.0
>>>> pull-group2-pbcatom     = 1343
>>>> pull_coord1_groups      = 1 2
>>>> pull-coord1-vec         = 0 0 1
>>>> pull_coord1_rate        = 0.0
>>>> pull_coord1_k           = 1000          ; kJ mol^-1 nm^-2
>>>> pull_coord1_start       = yes           ; define initial COM distance >
>>>> 0
>>>>
>>>> Segfaults without error messages within 20-50 ps. I suspected that
>>>> pull-coord1-vec is somehow hating zero coord1_rate, so I changed that to
>>>> 1.0e-9.
>>>>
>>>> The result is:
>>>>
>>>> ***
>>>> Program:     gmx mdrun, version 2018
>>>> Source file: src/gromacs/mdlib/sim_util.cpp (line 776)
>>>> MPI rank:    0 (out of 4)
>>>>
>>>> Fatal error:
>>>> Step 40700: The total potential energy is -nan, which is not finite.
>>>> The LJ
>>>> and electrostatic contributions to the energy are 15121 and -108386,
>>>> respectively. A non-finite potential energy can be caused by overlapping
>>>> interactions in bonded interactions or very large or Nan coordinate
>>>> values.
>>>> Usually this is caused by a badly- or non-equilibrated initial
>>>> configuration,
>>>> incorrect interactions or parameters in the topology.
>>>> ***
>>>>
>>>> No indication of exploding system at step 40000, all energies looking
>>>> legit.
>>>>
>>>> Any thoughts?
>>>>
>>>
>>> Looks buggy. Can you upload a .tpr to Redmine with a description of the
>>> problem?
>>>
>>> -Justin
>>>
>>> Alex
>>>>
>>>>
>>>>
>>>>
>>>> On Wed, Apr 18, 2018 at 4:00 PM, Alex <nedomacho at gmail.com> wrote:
>>>>
>>>> ... and everything is segfaulting, though worked perfectly fine until
>>>>> the
>>>>> 'cylinder' thing.
>>>>>
>>>>> On Wed, Apr 18, 2018 at 3:45 PM, Alex <nedomacho at gmail.com> wrote:
>>>>>
>>>>> Which then would be = 0 0 1, correct? This though will have to coexist
>>>>>> with pull_coord1_rate        = 0.0
>>>>>>
>>>>>> Is that reasonable?
>>>>>>
>>>>>> Thanks,
>>>>>>
>>>>>> Alex
>>>>>>
>>>>>> On Wed, Apr 18, 2018 at 2:11 PM, Justin Lemkul <jalemkul at vt.edu>
>>>>>> wrote:
>>>>>>
>>>>>>
>>>>>>> On 4/18/18 4:08 PM, Alex wrote:
>>>>>>>
>>>>>>> In addition, we seem to have a different problem:
>>>>>>>>
>>>>>>>> Cylinder geometry requires the use of pull-coord1-vec, not
>>>>>>> pull-coord1-dim.
>>>>>>>
>>>>>>> -Justin
>>>>>>>
>>>>>>>
>>>>>>> -------------------------------------------------------
>>>>>>>
>>>>>>>> Program:     gmx grompp, version 2018
>>>>>>>> Source file: src/gromacs/gmxpreprocess/readpull.cpp (line 234)
>>>>>>>>
>>>>>>>> Fatal error:
>>>>>>>> With pull geometry cylinder the pull vector can not be 0,0,0
>>>>>>>>
>>>>>>>>
>>>>>>>> On Wed, Apr 18, 2018 at 6:34 AM, Justin Lemkul <jalemkul at vt.edu>
>>>>>>>> wrote:
>>>>>>>>
>>>>>>>>
>>>>>>>> On 4/18/18 3:17 AM, Alex wrote:
>>>>>>>>>
>>>>>>>>> I suppose this question is mostly for Justin...
>>>>>>>>>
>>>>>>>>>> Let me remind what I am dealing with and ask if my idea is
>>>>>>>>>> correct.
>>>>>>>>>>
>>>>>>>>>> I have a rectangular membrane in XY with a pore at (X/2, Y/2) in
>>>>>>>>>> water
>>>>>>>>>> and want to get the Gibbs free energy curve for an ion. For this,
>>>>>>>>>> I
>>>>>>>>>> have a
>>>>>>>>>> bunch of starting configurations at (X/2, Y/2) and Z varying
>>>>>>>>>> between
>>>>>>>>>> some
>>>>>>>>>> -z0 and z0. The bias in the "fake" pull mdp is applied as N N Y.
>>>>>>>>>> Near
>>>>>>>>>> the
>>>>>>>>>> membrane, this means the entire plane is sampled, which adds
>>>>>>>>>> contributions
>>>>>>>>>> I am not interested in. I want the pore and a small region of the
>>>>>>>>>> membrane
>>>>>>>>>> around it, not the membrane, given its properties.
>>>>>>>>>>
>>>>>>>>>> So, I applied weak (k=50) in-plane restraint to the ion for each
>>>>>>>>>> of
>>>>>>>>>> the
>>>>>>>>>> sampled configurations -- to keep the sampling region a bit
>>>>>>>>>> closer to
>>>>>>>>>> the
>>>>>>>>>> pore. The results look completely different, but they finally make
>>>>>>>>>> very
>>>>>>>>>> good qualitative sense. The ion still walks around within a small
>>>>>>>>>> disk, but
>>>>>>>>>> not much -- this tentatively makes me happy.
>>>>>>>>>>
>>>>>>>>>> Would you believe the results obtained this way, assuming
>>>>>>>>>> otherwise
>>>>>>>>>> proper setup?
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> This is the exact situation that the "cylinder" pull geometry is
>>>>>>>>>>
>>>>>>>>> designed
>>>>>>>>> to model. I suggest going that route rather than ad hoc restraints,
>>>>>>>>> which
>>>>>>>>> will influence the PMF and should be otherwise accounted for.
>>>>>>>>>
>>>>>>>>> -Justin
>>>>>>>>>
>>>>>>>>> --
>>>>>>>>> ==================================================
>>>>>>>>>
>>>>>>>>> Justin A. Lemkul, Ph.D.
>>>>>>>>> Assistant Professor
>>>>>>>>> Virginia Tech Department of Biochemistry
>>>>>>>>>
>>>>>>>>> 303 Engel Hall
>>>>>>>>> 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
>>>>>>> Virginia Tech Department of Biochemistry
>>>>>>>
>>>>>>> 303 Engel Hall
>>>>>>> 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
> Virginia Tech Department of Biochemistry
>
> 303 Engel Hall
> 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.
>

More information about the gromacs.org_gmx-users mailing list