[gmx-users] Minimising forces for vibrational normal mode analysis

Peter Stern peter.stern at weizmann.ac.il
Tue Mar 1 04:16:29 CET 2016


This is a non-problem.  You will never get zero eigenvalues for such a large system since you are not at an absolute minimum.  The fact that you got no negative eigenvalues and the six lowest eigenvalues are very small (surely much, much smaller than the subsequent eigenvalues) indicates that you are at a reasonably good minimum. But don't be misled. Another (different) minimization could bring the system to a different minimum (also good) with somewhat different eigenvalues and eigenvectors.

Regards,
Peter Stern

Sent from my iPhone

> On 29 Feb 2016, at 3:23 PM, Nash, Anthony <a.nash at ucl.ac.uk> wrote:
> 
> Hi all,
> 
> Hoping this will be the final stumbling block. I’ve successfully minimised
> my structure to a a very small max force:
> Potential Energy  = -1.93736854460113e+03
> Maximum force     =  7.39399083846082e-04 on atom 34
> Norm of force     =  1.64068483698544e-04
> 
> 
> When I perform integrator=nm there are no warning of potential imaginary
> (negative) values. However, upon running:
> g_nmeig_d -f dogdic_nma.mtx -s dogdic_nma.tpr -last 3000
> I see the following:
> 
> 
> ------------------------------------
> Full matrix storage format, nrow=279, ncols=279
> Diagonalizing to find vectors 1 through 279...
> 
> One of the lowest 6 eigenvalues has a non-zero value.
> 
> This could mean that the reference structure was not
> properly energy minimized.
> Writing eigenvalues...
> 
> Back Off! I just backed up eigenval.xvg to ./#eigenval.xvg.2#
> Writing eigenfrequencies - negative eigenvalues will be set to zero.
> ------------------------------------
> 
> 
> 
> My lowest 6 eigenvalues are:
> 1       0.0780823
>     2        0.153357
>     3        0.174772
>     4        0.327362
>     5        0.466203
>     6         0.64693
> 
> 
> I’m a little confused by the statement "One of the lowest 6 eigenvalues
> has a non-zero value. This could mean that the reference structure was not
> properly energy minimized.” Really? Surely it out to be if one or more of
> the values is negative.  Also, there were no eigenvalues set to zero
> (hence, I can only assume I have no negative eigenvalues).
> 
> Would appreciate a little insight.
> Many thanks
> Anthony
> 
> Dr Anthony Nash
> Department of Chemistry
> University College London
> 
> 
> 
> 
> 
> On 29/02/2016 16:25, "gromacs.org_gmx-users-bounces at maillist.sys.kth.se on
> behalf of Nash, Anthony"
> <gromacs.org_gmx-users-bounces at maillist.sys.kth.se on behalf of
> a.nash at ucl.ac.uk> wrote:
> 
>> Dear Mark and Justin,
>> 
>> By removing the restraints (your suggestions) it appears to have worked!
>> 
>> Maximum force: 9.21098e-04
>> Writing Hessian...
>> 
>> This now matches the same force as my energy minimisation. Many thanks for
>> you help.
>> 
>> Thanks
>> Anthony
>> 
>> 
>> 
>> Dr Anthony Nash
>> Department of Chemistry
>> University College London
>> 
>> 
>> 
>> 
>> 
>> On 29/02/2016 14:46, "gromacs.org_gmx-users-bounces at maillist.sys.kth.se on
>> behalf of Mark Abraham" <gromacs.org_gmx-users-bounces at maillist.sys.kth.se
>> on behalf of mark.j.abraham at gmail.com> wrote:
>> 
>>> Hi,
>>> 
>>> An earlier mdp file suggested you were using position restraints. There
>>> should be no need for this, nor any problem, but what happens without
>>> them?
>>> 
>>> Mark
>>> 
>>>> On Mon, 29 Feb 2016 15:39 Nash, Anthony <a.nash at ucl.ac.uk> wrote:
>>>> 
>>>> Hi Justin,
>>>> 
>>>> After some digging I had found that link and made some adjustments (as
>>>> presented in the later email).
>>>> 
>>>> After a series of energy minimisations (including switching LINCS off,
>>>> and
>>>> dropping the energy step to a very small number), and with the final
>>>> command:
>>>> grompp_d -f cg.mdp -c modic_cg_3.gro -p system.top -o modic_cg_4 -t
>>>> modic_cg_3.trr
>>>> 
>>>> 
>>>> Polak-Ribiere Conjugate Gradients converged to Fmax < 0.001 in 6883
>>>> steps
>>>> Potential Energy  = -1.92681826422996e+03
>>>> Maximum force     =  8.54361662283108e-04 on atom 44
>>>> Norm of force     =  2.71110116926712e-04
>>>> 
>>>> 
>>>> However, when switching to integrator=nm, running grompp_d -f nma.mdp
>>>> -c
>>>> modic_cg_4.gro -p system.top -o modic_nma -t modic_cg_4.trr, and then
>>>> mdrun_d, I get:
>>>> 
>>>> Maximum force: 5.19179e+01
>>>> The force is probably not small enough to ensure that you are at a
>>>> minimum.
>>>> Be aware that negative eigenvalues may occur
>>>> when the resulting matrix is diagonalized.
>>>> 
>>>> 
>>>> 
>>>> I¹m still struggling to yield a maximum force during the integrator=nm
>>>> step, as presented from the earlier cg step. The .mdp files are
>>>> identical
>>>> with the exception of the integrator line.
>>>> 
>>>> Thanks
>>>> Anthony
>>>> 
>>>> 
>>>> Dr Anthony Nash
>>>> Department of Chemistry
>>>> University College London
>>>> 
>>>> 
>>>> 
>>>> 
>>>> 
>>>> On 29/02/2016 12:49, "gromacs.org_gmx-users-bounces at maillist.sys.kth.se
>>>> on
>>>> behalf of Justin Lemkul"
>>>> <gromacs.org_gmx-users-bounces at maillist.sys.kth.se on behalf of
>>>> jalemkul at vt.edu> wrote:
>>>> 
>>>>> 
>>>>> 
>>>>>> On 2/29/16 3:41 AM, Nash, Anthony wrote:
>>>>>> Hi Tsjerk,
>>>>>> 
>>>>>> The two .mdp files are virtually identical (the only exception being
>>>>>> what
>>>>>> defines one as a conjugate-gradient, and the other for normal mode
>>>>>> analysis):
>>>>>> 
>>>>>> CONJUGATE GRADIENT:
>>>>>> define = -DPOSRES
>>>>>> integrator      = cg
>>>>>> emtol           = 0.001
>>>>>> emstep          = 0.0002
>>>>>> nsteps          = 1000000
>>>>>> nstcgsteep      = 100
>>>>>> cutoff-scheme = verlet
>>>>>> nstlist         = 10
>>>>>> ns_type         = grid
>>>>>> rlist           = 1.4
>>>>>> coulombtype     = PME
>>>>>> rcoulomb        = 1.4
>>>>>> rvdw            = 1.4
>>>>>> pbc             = xyz
>>>>>> 
>>>>>> 
>>>>>> NORMAL MODE ANALYSIS
>>>>>> define = -DPOSRES
>>>>>> integrator = nm
>>>>>> emtol = 0.001
>>>>>> emstep = 0.0002
>>>>>> nsteps = 1000000
>>>>>> cutoff-scheme = verlet
>>>>>> nstlist = 10
>>>>>> ns_type = grid
>>>>>> rlist = 1.4
>>>>>> coulombtype = PME
>>>>>> rcoulomb = 1.4
>>>>>> rvdw = 1.4
>>>>>> pbc = xyz
>>>>>> 
>>>>>> 
>>>>>> The cg energy minimisation did NOT result in any warning about force
>>>> not
>>>>>> converging. The result I got (I just re did it now) was:
>>>>>> Polak-Ribiere Conjugate Gradients converged to Fmax < 0.001 in 6994
>>>>>> steps
>>>>>> Potential Energy  = -1.73087278108256e+03
>>>>>> Maximum force     =  9.97199385029354e-04 on atom 72
>>>>>> Norm of force     =  4.63373534634654e-04
>>>>>> 
>>>>>> 
>>>>>> But then when I run normal mode analysis (integrator=nm) I get:
>>>>>> Maximum force: 6.97334e+02
>>>>>> The force is probably not small enough to ensure that you are at a
>>>>>> minimum.
>>>>>> Be aware that negative eigenvalues may occur
>>>>>> when the resulting matrix is diagonalized.
>>>>>> 
>>>>>> 
>>>>>> My work flow:
>>>>>> grompp_d -f cg.mdp -c modic_en.gro -p system.top -o modic_cg
>>>>>> mdrun_d -deffnm modic_cg
>>>>>> grompp_d -f nma.mdp -c modic_cg.gro -p system.top -o modic_nma
>>>>>> mdrun_d -deffnm modic_nma
>>>>> 
>>>>> Relying only on .gro format is insufficient precision to continue from
>>>>> minimization to NMA.  Make sure to pass the .trr to grompp -t.
>>>>> 
>>>>> http://www.gromacs.org/Documentation/How-tos/Normal_Mode_Analysis
>>>>> 
>>>>> -Justin
>>>>> 
>>>>> --
>>>>> ==================================================
>>>>> 
>>>>> Justin A. Lemkul, Ph.D.
>>>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>>>> 
>>>>> Department of Pharmaceutical Sciences
>>>>> School of Pharmacy
>>>>> Health Sciences Facility II, Room 629
>>>>> University of Maryland, Baltimore
>>>>> 20 Penn St.
>>>>> Baltimore, MD 21201
>>>>> 
>>>>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>>>>> http://mackerell.umaryland.edu/~jalemkul
>>>>> 
>>>>> ==================================================
>>>>> --
>>>>> 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.
>>>> 
>>>> --
>>>> 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.
>>> -- 
>>> 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.
>> 
>> -- 
>> 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.
> 
> -- 
> 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