[gmx-users] hints for core/shell optimization?

[Justin Lemkul]

On 8/1/14, 6:13 AM, Tamas Karpati wrote:
> Dear Justin,
>
>>> If I have a Mg++ ion with q_D=+3.2 then for the shell particle I'd set
>>> +3.2 and for the Mg atom I should stay with +2? So far -in order to maintain
>>> charge neutrality of the system- I've been using -1.2 for Mg to have the
>>> correct +2e for the AD duett. Big difference!
>>
>> Yes, the Drude and core atom will have charges of opposite sign to yield the
>> total charge on the pair.  A positively charged auxiliary particle is a bit
>> weird, but not wrong, if the model was parametrized that way.  Normally the
>
> I've tried writing +2 charge for Mg in the TOP file and let GROMACS combine
> it with the Drude charge of +3.2 but it has set a total charge of ca.
> few thousands
> e and induced a great movement of the shell particles. The average A-D
> distance raised from ca. 0.01 to 0.3 nm! I consider the original solution
> better (I pre-combine q_D with q_A and write q_C=q_A-q_D into the TOP file).
>

Not sure I follow (excerpts from the .top are preferred); my suggestion was not 
to assign +2 to the Mg2+ core atom, as that would not make sense.  The value of 
q_total (q_A + q_D) should be 2, so the original approach was right.  I was only 
questioning the logic of having a positively charged auxiliary particle.  That 
is quite exotic (though, again, not necessarily wrong).

>> The way I have done the exclusions is redundant and in principle not totally
>> necessary.  Atom-Drude pairs connected via bonds or within [polarization]
>> directives (e.g. 1-2 and 3-4) should already be excluded from one another,
>> so in reality:
>>
>> [ exclusions ]
>> 1 3 4
>> 2 3 4
>>
>> should cover it.
>
> I does, at least my job is functionally runs but MDRUN stops for too low
> atomic displacements while forces are in the sky:
>
> # Step=    2, Dmax= 5.0e-03 nm, Epot= -1.35892e+08 Fmax= 3.49164e+07, atom= 1242
> # Step=    3, Dmax= 6.0e-03 nm, Epot= -1.72899e+08 Fmax= 8.98685e+08, atom= 1900
> # Step=    7, Dmax= 9.0e-04 nm, Epot= -3.21500e+09 Fmax= 1.05629e+16, atom= 6794
> # Step=   18, Dmax= 1.1e-06 nm, Epot= -1.33283e+09 Fmax= 1.08140e+14, atom= 679
> # Energy minimization has stopped, but the forces have not converged to the
> # requested precision Fmax < 0.1 (which may not be possible for your system).
> # It stopped because the algorithm tried to make a new step whose size was too
> # small, or there was no change in the energy since last step. Either way, we
> # regard the minimization as converged to within the available machine
> # precision, given your starting configuration and EM parameters.
>
> 1. Usual landing at ca. 1e-6 nm for Dmax...OK
> 2. Last Epot is usually higher than the one before...?
> 3. Fmax almost monotonically increases...?
> 4. Not all steps are written, although I'd like to look at them...?

This means that intermediate steps were unsuccessful; mdrun only writes an 
update when the potential energy has decreased, so for instance, steps 9-17 
produced no useful modification.

> (After MDRUN -using dmxdump and editconf- I put normal atoms
> and shell particles in a separate animated XYZ files to see morphing
> the input into an output and somtimes miss more details. If it was possible,
> I'd take a look also at the Drude optimization steps between two EM steps.)
>
> Setting a_Thole to a different value or rescale too large literature
> polarizabilities
> do not help. In fact, there is no big difference even between applying
> or omitting
> [thole_polarization] along with [polarization].
>

Can you provide a link to the paper(s) that describe your model?  I've been 
working largely blind here trying to guess at what the issues might be, but if 
Thole screening isn't making a big impact, then it may not be part of the model 
at all.

> One more -for me, now unfortunate- consequence of Thole-ing the
> polarizability: after successfully avoiding through-cell bonds in my
> PBC simulation, I've now re-introduced them via screening.
> Had to set "periodic-molecules = yes" in the MDP file, yet have
> the following question.
>
> There are big many A-D particle pairs very close to the walls and,
> of course, to other dipoles on the other side.
> What if I simply Thole+exclude them like fully internal AD_i-AD_j
> pairs of particle pairs? Will it be a long quasi-bond from one wall
> to the opposide side-wall or it will just be a short through-wall bond?
>

If you are using periodic_molecules = yes, then it will be a bond through the 
wall, not across the unit cell.

> I think the majority of the pairs of AD combos are fully inside the
> box so I don't expect the full solution from this but it needs be clarified.
>
> You've mentioned correcting the Thole-code by a scale factor.
> Can I emulate this correction by scaling the polarizabilities
> in the TOP file?

No, I've made wholesale changes to the code in the way that energies and forces 
are computed.  The current code has bugs.  If you want a modified version (beta 
software, so caveat emptor), contact me off-list.

-Justin

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

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

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