Re: [gmx-users] hints for core/shell optimization? (Tamas Karpati)
Dear Tamas, Just a small note on the subject. Keep in mind that particular choice of atom and shell values (while q_A+q_D remains constant) may influence some properties of the system. More details can be found in Asmadi at al. (Journal of Molecular Liquids 188 (2013) 245-251). Although, for ions this effect, most probably, is not very large, it is worth to take it into account. Best regards, Mikhail Mikhail Stukan, Ph.D. Senior Research Scientist Schlumberger Dhahran Carbonate Research Center Dhahran Techno Valley - KFUPM P.O. Box 39011, Dammam/Doha Camp 31942 Kingdom of Saudi Arabia Tel: +966 3 331 6182 Fax: +966 3 330 0845 mstu...@slb.commailto:mstu...@slb.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-requ...@gromacs.org.
Re: [gmx-users] hints for core/shell optimization?
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
Re: [gmx-users] hints for core/shell optimization?
Dear Justin, I'm in the experimenting with [exclusions] for [thole_polarization]. Still not there but approaching... I'd like to ask more help as the Manual is not detailed enough for me to comprehend the whole. When using the core/shell model I expect the following extra effects and functions beyond the normal simulations. [Denotations: A = core atom, D = Drude (or shell) particle, AD_i = the i_th A-D oscillator (dipole).] 1. With a q_D charge on D, the normal q_A becomes q_C=q_A-q_D to preserve the q_A=q_C+q_D atomic charge in a c/s simulation (since D is attached to A). I assume that GROMPP expects q_C for atoms, q_D for shell particles in the GRO files. Is this correct? 2. I assume that GROMACS evaluates _Coulomb_potentials_ only between atoms and only using q_A values, while q_D never ever occurs in such expressions. (Of course, technically q_A is obtained as q_C+q_D.) Is this correct or -on top of each (Ai,Aj) pair- all (Di,Dj), (Di,Aj) and (Ai,Dj) contributions are added to the above Coulomb sum? (See also the merge-comment below 3B.) 3. I expect that the _polarization_potential_ is calculated by either one of the following 2 possibilities (A or B) after the D positions have been optimized in each EM or MD step (I guess forces are evaluated similarly for that optimization). A) The AD_i - AD_j pair potentials are summed using q_D and -q_D values for the i_th and j_th dipoles. For a given pair one has Coulombic terms for the (Ai,Aj), (Ai,Dj), (Aj,Di) and the (Di,Dj) charged point pairs. I assume that the (Ai,Di) and (Aj,Dj) terms are only used in the shell particle optimizations. [Nothing is added twice, eg. in (Ai,Aj), as the normal (non c/s) part used the q_A charges for the Coulomb interaction, while here we have +/-q_D charges.] B) The mu_i and mu_j dipole moments of the AD_i - AD_j pairs are summed (of course, intrinsically calculated from the i_th and j_th q_D charges and A-D distances). Is A or B correct? Which one? (I can imagine this polarization part merged with the Coulomb part in 2. for efficiency but I'd like to make the conceptual distinction.) 4. If the i_th and j_th dipoles are too close to one another then the polarization becomes catastrophic which is indicated by too high Coulombic forces, mu_AD dipole moments or A-D distances (not infinite polarizabilities as alpha is an input variable). Here is the entry point for Thole screening: for dipoles in one another's proximity scale down the Coulomb forces between just the two Drude (shell) particles. The scale factor is 1-(1+r12/2)*exp(-r12), where r12=a*d/(alpha_1*alpha_2)^1/6 with a: magic number, calibrated for each material, d: distance of D_i and D_j and alpha_x: polarizabilities [nm^3]. Thus shell-shell (Di,Dj) catastrophe have been avoided but what is with the remaining (Ai,Aj), (Ai,Dj), (Aj,Di) terms? (I have tried using the [excludes] directive to remove these but my simulation still couldn't converge. I also excluded the Di-Dj, Ai-Di and Aj-Dj pairs with not more success. I'm not even convinced that the inner-loop optimization for the shell particles completes so I've set niter=1.) Also confusing is that -looking at the MDRUN output and the Fmax values that increase from the initial 1e4 to even 1e16 during the EM process- I notice both normal atoms and Drude particles having the highest force acting on them in one or the other EM step. Besides, Epot is gradually decreasing so I'm getting more happy with the results. However, full convergence is not achieved. My goal is to learn what to [exclude] and [thole_polarization] in my TOP file after locating pairs of AD particle doublets that are too close to each other. I imagine that for each AD_i-AD_j pairs I need to do the same. Could you please shed more light on what same means here? I thank you for your continuous support. With regards, toma PS: invented a work-around for GROMPP to accept the [exclusions] directive in the TOP file when no bonds are present (the crystal is held together by Buckingham+Coulomb): one needs to just write an empty [bonds] directive (without a single bond specified). Otherwise it claims the [exclusions] is at a wrong position in TOP. -- 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-requ...@gromacs.org.
Re: [gmx-users] hints for core/shell optimization?
Dear Justin, 1. With a q_D charge on D, the normal q_A becomes q_C=q_A-q_D to preserve the q_A=q_C+q_D atomic charge in a c/s simulation (since D is attached to A). I assume that GROMPP expects q_C for atoms, q_D for shell particles in the GRO files. Is this correct? No. Every atom in the topology is assigned its charge. grompp doesn't manipulate them. You need to assign q_A and q_D, e.g. as is done in sw.itp. The contents of [polarization] are only used to calculate the force constant of the bond between the atom and Drude (see equation in sw.itp or code in src/gromacs/gmxlib/bondfree.c) 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! Anything with a charge contributes to the Coulomb potential. OK, I see now why to [exclude] lots of things. Neither. The polarization energy is the output of the polarize() function (again, see src/gromacs/bondfree.c) and is a simple bonded potential using the atom - Drude distance and the value of k (force constant) that is computed from q_D and alpha. Thanks. I'm going to check the source code. Thole screening is applied to all possible pairs. Remember, it's With possible meaning close-enough dipoles? dipole-dipole screening, so all of those terms are screened, hence the I think I took words too sharply... I begin to realize that dip-dip screening means point charge-point charge screening. Thus, eg. with AD_i and AD_j close enough we have 4 particles and 4*3/2=6 point-to-point interactions to screen. (I did it, not yet there.) reason you need to specify 4 atoms in the [thole_polarization] directive. See thole_pol() and do_1_thole() in bondfree.c. For once, the CHARMM documentation is actually better than Gromacs (though I am in the process of fixing that): http://www.charmm.org/documentation/c34b1/drude.html I'm sure it worths looking at. Appreciated. This still suggests topological instability. Or, as I said before, incorrect force calculation that is in do_1_thole() that I found to cause failures in otherwise sensible input files for our FF. I'm looking forward to seeing the update :) Exclude nearest neighbors, apply Thole screening to them. Depending on the geometry of the system, perhaps second neighbors, as well. For instance, if atom 1 is an atom, 2 is its Drude, 3 is an atom, and 4 is its Drude: [ exclusions ] 1 2 3 4 2 1 3 4 3 1 2 4 4 1 2 3 [ thole_polarization ] ; ai aj ak al func etc. 1 2 3 42 a alpha_12 alpha_34 Ooops... I expected finding something like that in a tutorial. But I see that you specify pairs twice: 1-2 and 2-1 etc. Maybe I just eliminated half of the polarization and screening potentials? With regards, toma -- 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-requ...@gromacs.org.
Re: [gmx-users] hints for core/shell optimization?
On 7/31/14, 2:17 PM, Tamas Karpati wrote: Dear Justin, 1. With a q_D charge on D, the normal q_A becomes q_C=q_A-q_D to preserve the q_A=q_C+q_D atomic charge in a c/s simulation (since D is attached to A). I assume that GROMPP expects q_C for atoms, q_D for shell particles in the GRO files. Is this correct? No. Every atom in the topology is assigned its charge. grompp doesn't manipulate them. You need to assign q_A and q_D, e.g. as is done in sw.itp. The contents of [polarization] are only used to calculate the force constant of the bond between the atom and Drude (see equation in sw.itp or code in src/gromacs/gmxlib/bondfree.c) 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 auxiliary particle represents deformation of the electron cloud, so it is negative while the core atom is positive. Anything with a charge contributes to the Coulomb potential. OK, I see now why to [exclude] lots of things. Neither. The polarization energy is the output of the polarize() function (again, see src/gromacs/bondfree.c) and is a simple bonded potential using the atom - Drude distance and the value of k (force constant) that is computed from q_D and alpha. Thanks. I'm going to check the source code. Small correction (sorry for the typo) - the file is src/gromacs/gmxlib/bondfree.c. Thole screening is applied to all possible pairs. Remember, it's With possible meaning close-enough dipoles? dipole-dipole screening, so all of those terms are screened, hence the I think I took words too sharply... I begin to realize that dip-dip screening means point charge-point charge screening. Thus, eg. with AD_i and AD_j close enough we have 4 particles and 4*3/2=6 point-to-point interactions to screen. (I did it, not yet there.) reason you need to specify 4 atoms in the [thole_polarization] directive. See thole_pol() and do_1_thole() in bondfree.c. For once, the CHARMM documentation is actually better than Gromacs (though I am in the process of fixing that): http://www.charmm.org/documentation/c34b1/drude.html I'm sure it worths looking at. Appreciated. This still suggests topological instability. Or, as I said before, incorrect force calculation that is in do_1_thole() that I found to cause failures in otherwise sensible input files for our FF. I'm looking forward to seeing the update :) It will be a while. It works, but is not yet compatible with DD (coding ongoing) and needs to be tested quite a bit. Hopefully soon... Exclude nearest neighbors, apply Thole screening to them. Depending on the geometry of the system, perhaps second neighbors, as well. For instance, if atom 1 is an atom, 2 is its Drude, 3 is an atom, and 4 is its Drude: [ exclusions ] 1 2 3 4 2 1 3 4 3 1 2 4 4 1 2 3 [ thole_polarization ] ; ai aj ak al func etc. 1 2 3 42 a alpha_12 alpha_34 Ooops... I expected finding something like that in a tutorial. But I see that you specify pairs twice: 1-2 and 2-1 etc. Maybe I just eliminated half of the polarization and screening potentials? 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. -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 jalem...@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-requ...@gromacs.org.
Re: [gmx-users] hints for core/shell optimization?
Dear Justin, No, they're not mutually exclusive. The Thole function screens nearby dipoles so they don't experience artificial forces. The thole_polarization name is a bit unfortunate; it should really be thole_screening instead. I was probably confused by the two contribs in the MDRUN log file. Thank you for the clarification. I suspect you're seeing large forces due to one (or more effects). If your model doesn't have any bonds, then you don't have any excluded interactions, so instead of getting a screened Coulombic interaction via Thole screening, you're getting both added together. You should make use of [exclusions] for neighboring pairs, likely, but I guess the bigger question is how your force field was parametrized and what its expected behavior should be. Do you mean that for the Thole screening to work I need bonds in order to have 1-2- and 1-3 interactions? Or do I need the exclusions that pull bonds in? Should I declare bonds of an average (eg. 0.2 nm) length with k=0 force constants then add [exclusions] for farther atoms? Please confirm what I need to exclude (the neighboring, 1-2, pairs rather than the 1-3 or farther pairs?) FYI: I took the model and parameters (fitted to experimental geometries and diffusion enthalpies) from a (non-GROMACS) paper on some crystal. The model seems chaotic to me but is perfect to learn MM -coming from the QC world. The other effect is one that I noticed this weekend. I think the forces coming out of the Thole function are incorrect; at least they are in the case of our Drude polarizable force field. I changed the force calculation and everything works as expected now, in my development version of the code. Nice to hear that. I'm looking forward to seeing it at your site for download :) field parameters? In our Drude model, 1-2 and 1-3 interactions (based on atom connectivity, not Drude connectivity) are excluded and electrostatics calculated via Thole. You don't have any bonds, so the choice here is not exactly clear and I don't want to pre-judge without knowing what the force field is or how it should be expected to behave. As I said my FF is built from literature and I consider it chaotic. The crystal is held together by nonbonding interactions and polarization forces. In addition, I do also have some angle forces (without any bonds). The core/shell model means just some extra shell particles on selected atom types; c/s -s are bonded by [polarization] items. Recently -with your aid- I started screening them by the Thole method. (Of course, I did otherwise but I will correct my method, thanks :) No dummies are applied. Although I though about it, I haven't found in the Manual a way to put D on a single A (rather than on 2 or more atoms). In such a way I could have reduced charges on the atoms so that q_D = -q_S and [polarization] would have D-S instead of A-S bonds. Only I couldn't figure out how to fix D on A sites. Is Drude connectivity something I also use (no bonds)? In a sense I then have both Drude and Buckingham connectivities. Since undocumented, I do not really know in which (atom-atom, atom-shell and shell-shell) cases do we have Coulombic, point charge-dipole and dipole-dipole interactions and whether we can switch them on/off. Finally, only I expect from the c/s model to make atom-atom interactions more sophisticated: they have an extra channel of communication -between their polarizable aspect. GROMACS will crash if I put shells exactly at the atomic positions thus I use some 0.01 nm offset (in random direction to avoid metastable symmetric structures). I do not expect catastrophic polarization with such small offsets (small compared to typical bond lengths of 0.2 nm) but something goes wrong as both shell and atomic optimizations (EM run) stop after a few steps without approaching the force criterium. As it should; the value of a should be the sum of the two atomic Thole factors, or 2.6, whichever your force field uses. Thole screening was not mentioned in the paper I started from. I started from a=2.6 which was suggested in the GROMACS Manual. Later I noticed that this value was estabilished for ethanol and such. I've seen values between ca. 0.5 and 2.6 so I played between 0.1 and 10 to see that it has a strong effect on forces and the final structure and see that it is to be calibrated for a new material. I understand that a controls the drop of interaction with distance. Maybe I need to use unusual values. Sounds like it could be incorrect forces, but without seeing the Thole section of your topology, I can't say for sure. The Thole section in my TOP file lists all c/s pairs that are closer than R (eg. 0.5 nm) to one another (also specifying func=2, a and alpha values). 0.5 - 1.0 nm sounds like far too large of a radius for Thole interactions. Should I use something like a range [0.3;0.5]? (Typical bonds are shorter than 0.3 nm.) Or should I use [0;0.3] to just save the neighbors from the
Re: [gmx-users] hints for core/shell optimization?
On 7/30/14, 5:14 AM, Tamas Karpati wrote: Dear Justin, No, they're not mutually exclusive. The Thole function screens nearby dipoles so they don't experience artificial forces. The thole_polarization name is a bit unfortunate; it should really be thole_screening instead. I was probably confused by the two contribs in the MDRUN log file. Thank you for the clarification. I suspect you're seeing large forces due to one (or more effects). If your model doesn't have any bonds, then you don't have any excluded interactions, so instead of getting a screened Coulombic interaction via Thole screening, you're getting both added together. You should make use of [exclusions] for neighboring pairs, likely, but I guess the bigger question is how your force field was parametrized and what its expected behavior should be. Do you mean that for the Thole screening to work I need bonds in order to have 1-2- and 1-3 interactions? Or do I need the exclusions that pull bonds in? Should I declare bonds of an average (eg. 0.2 nm) length with k=0 force constants then add [exclusions] for farther atoms? Please confirm what I need to exclude (the neighboring, 1-2, pairs rather than the 1-3 or farther pairs?) FYI: I took the model and parameters (fitted to experimental geometries and diffusion enthalpies) from a (non-GROMACS) paper on some crystal. The model seems chaotic to me but is perfect to learn MM -coming from the QC world. I'm hesitant to start making suggestions about modifications without knowing the exact details of the model. There are often aspects left unexplained in such papers, and I suspect there are special aspects of whatever software the authors used that control these physical properties. If they didn't use any screening function explicitly, then their code must be already handling it somehow. But in short, no, bonds are not compulsory for using Thole, but exclusions definitely are. I would exclude at least the nearest-neighbor interactions. Can you post a snippet of your [polarization] and [thole_polarization] directives? The other effect is one that I noticed this weekend. I think the forces coming out of the Thole function are incorrect; at least they are in the case of our Drude polarizable force field. I changed the force calculation and everything works as expected now, in my development version of the code. Nice to hear that. I'm looking forward to seeing it at your site for download :) field parameters? In our Drude model, 1-2 and 1-3 interactions (based on atom connectivity, not Drude connectivity) are excluded and electrostatics calculated via Thole. You don't have any bonds, so the choice here is not exactly clear and I don't want to pre-judge without knowing what the force field is or how it should be expected to behave. As I said my FF is built from literature and I consider it chaotic. The crystal is held together by nonbonding interactions and polarization forces. In addition, I do also have some angle forces (without any bonds). The core/shell model means just some extra shell particles on selected atom types; c/s -s are bonded by [polarization] items. Recently -with your aid- I started screening them by the Thole method. (Of course, I did otherwise but I will correct my method, thanks :) No dummies are applied. Although I though about it, I haven't found in the Manual a way to put D on a single A (rather than on 2 or more atoms). In such a way I could have reduced charges on the atoms so that q_D = -q_S and [polarization] would have D-S instead of A-S bonds. Only I couldn't figure out how to fix D on A sites. Just a note for clarity here, since it might be a bit confusing. A Drude and a shell are the same thing; dummies (virtual sites) haven't been mentioned yet and don't sound like they're relevant. A-S is an atom-shell(Drude) bond and D-S is a dummy(virtual site)-shell bond. The use of S and D to mean different things when they could be the same in other nomenclature systems is somewhat confusing. Is Drude connectivity something I also use (no bonds)? In a sense I then have both Drude and Buckingham connectivities. What I meant was this (hopefully this shows up correctly) D1 D2 D3 ||| A1---A2---A3 The exclusions and Thole screening are based on atom connectivity only. Drudes D1 and D3 are technically separated by 4 bonds, but atoms A2 and A3 are only separated by 2, so the 1-3 interactions between A1, D1, A3, and D3 are all excluded and treated by the screening function. That is, a bond to a Drude/shell does not count as a bond when it comes to determining whether or not interactions should be excluded. Since undocumented, I do not really know in which (atom-atom, atom-shell and shell-shell) cases do we have Coulombic, point charge-dipole and dipole-dipole interactions and whether we can switch them on/off. Finally, only I expect from the c/s model to make atom-atom interactions more sophisticated: they
Re: [gmx-users] hints for core/shell optimization?
On 7/30/14, 9:37 AM, Tamas Karpati wrote: But in short, no, bonds are not compulsory for using Thole, but exclusions definitely are. I would exclude at least the nearest-neighbor interactions. Does this mean that by [thole_polarization] I define bonded interactions which can be dropped by [exclusions]? Would it not be easier to not define Thole at all for the nearest-neighbors in the first place? No, exclusions work on nonbonded interactions. Thole screening is, in theory, a nonbonded interaction, but it is handled as a bonded interaction in the code, which allows it to act on otherwise excluded atom pairs. Can you post a snippet of your [polarization] and [thole_polarization] directives? Sure. Shell particles directly follow their corresponding core (part. type A). Atoms 1 and 3 are metallic, 11 and 13 are nonmetallic elements. For the rest of the atom types (ions) there are no polarizations prescribed. [ polarization ] ; aiaj funct alpha 1 2 1 6.7 3 4 1 6.7 1112 1 1.5 1314 1 1.5 ... [ thole_polarization ] ; atom_i shell_i atom_j shell_j func a alpha_i alpha_j 1 2 3 4 2 2.6 6.7 6.7 1 2 11 12 2 2.6 6.7 1.5 ... I think your units are wrong here. The alpha values should be in nm^3. Those values of alpha are about 1000x larger than I think they should be. I assume it is irrelevant whether I use S or D particle type, thus I kept S. In my TOP file there are [moleculetype], [atoms], [angles], [polarization], [thole_polarization], [system] and [molecules] sections, nothing else. Keep S. A dummy (D) particle is a virtual site; very different! I quasi-manually put polarization on 2 types and put Thole on all pairs of A-S pairs that are within 0.35 nm. Should I then purge via [exclusions] all that are within one bond distance? Or should I just not put Thole there? Anything that is treated by Thole should be listed as an exclusion. The Manual sais The exclusions for non-bonded interactions are generated by grompp for neighboring atoms up to a certain number of bonds away... which is embarassing for two reasons: 1. [thole_polarization] is considered as bonded interaction I don't see how that is inconsistent with the manual. 2. I create GRO and TOP files then feed them (with MDP) to GROMPP to generate the TPR file. This is to control (=understand) every possible details of what's goind on. Do I miss something happening in the TOP --[GROMPP]-- TPR conversion? Just a note for clarity here, since it might be a bit confusing. A Drude and a shell are the same thing; dummies (virtual sites) haven't been mentioned yet and don't sound like they're relevant. A-S is an atom-shell(Drude) bond and D-S is a dummy(virtual site)-shell bond. The use of S and D to mean different things when they could be the same in other nomenclature systems is somewhat confusing. Thanks for clarifying that. What I meant was this (hopefully this shows up correctly) D1 D2 D3 ||| A1---A2---A3 The exclusions and Thole screening are based on atom connectivity only. Drudes D1 and D3 are technically separated by 4 bonds, but atoms A2 and A3 are only separated by 2, so the 1-3 interactions between A1, D1, A3, and D3 are all excluded and treated by the screening function. That is, a bond to a Drude/shell does not count as a bond when it comes to determining whether or not interactions should be excluded. Thanks for the figure. (Do you mean A1 and A3 are separated by 2 bonds?) Yes. As I understand, Thole scr. isn't based on but rather it creates a kind of atom connectivity which [exclusions] can partially remove. Eg. A1-D3 and D1-A3 get deleted for being 1-3 interactions. This is done by GROMPP automagically applying [exclusions] to my TOP file which doesn't have any. Have you verified this by using gmxdump? Exclusions are generated from bonds. You said you don't have any bonds, so there shouldn't be any exclusions aside from those you are adding. This, however, remained hidden from me (only appeared in the binary TPR file). That's what gmxdump is for. Then MDRUN's EM job applied [thole_polarization] rather than [polarization] on the excluded connections. On the other hand, the Manual explains in 4.4.3. Thole polarization that V_Thole is defined between 2 shell particles. In this case, I don't understand why to include atom indices in the [thole_polarization] section. I could simply define [polarization] for all my A-S pairs and [thole_polarization] for all the S-S pairs that are too close -and done. Thole screening is applied to atom-shell(Drude) pairs to account for unphysical dipole-dipole interactions. The manual is in need of updating. Even if I misunderstand something, for my simple scheme XYZ --[my code]- GRO + TOP --[GROMPP]-- TPR --[MDRUN]-- ... the automatisms in GROMACS seem to be overcomplicated. Thank you for helping me grabbing concepts and usage together:)
Re: [gmx-users] hints for core/shell optimization?
Dear Justin, I've tried applying the [thole_polarization] section in the format you've suggested. I've noticed the followings. 1. First I made a mistake and forgot to remove the [polarization] section from the TOP file. Even this way MDRUN worked and indicated the usual ca. 10 kJ/mol contribution from polarization. Besides a new contribution of ca. 600 kJ/mol appeared from Thole pol. I switched the former off to have just one kind of pol effect. Is this correct? 2. I assumed that in the [thole_polarization] section of the TOP file I needed to select which dipole pairs are to be included on basis of my chemical intuition -if I had any. I tried a few radii of sphere-of-inclusion but saw little differences. Is it a correct selection for Thole screening? (I have a crystal with separate ions held together by Buckingham forces and Thole polarization -there are no bonds at all.) Do I have a way to switch on/off intercell dipole interactions? 3. Playing with the a factor of the Thole scheme resulted in more pronounced changes in the contribution from polarization but I still see too strong Coulombic interactions. Omitting shell particles and thus polarization do scale Coulombics down by a factor of 100 or 1000. I guess that using [polarization] all dipoles interact with all the others (although no idea how dipoles interact through walls of the cell). I imagine that in the [thole_polarization] case only those dipole pairs will interact that I have mindlessly selected. In contrast with my expectation I noticed an order of magnitude greater Thole pol. values than Polarization values in the simple pol. case. (For your information, I have ca. 2000 atoms + 1000 shell particles in a box of ca. 3x3x3 nm. In case of the Thole scheme I have ca. 1e5 pcs. of dipole pairs with a 1 nm radius for inclusion while ca. 15000 for a radius of 0.5 nm.) The total potential is almost a pure Coulombic one and is too large. Do you have any hint on what could go wrong? Thanks for your attention. With best regards, toma On Fri, Jul 25, 2014 at 9:31 PM, Tamas Karpati tkarp...@gmail.com wrote: Dear Justin, Thank you for the valuable information. I'm starting experimenting with it. Have a nice weekend, toma On Fri, Jul 25, 2014 at 7:53 PM, Justin Lemkul jalem...@vt.edu wrote: On 7/25/14, 12:40 PM, Tamas Karpati wrote: Dear Justin, Thanks for your educational answers. Coulombic interactions fail at short distances; you probably need to apply I was afraid of that... somehow removing shells from the cores in the initial structure have let it functionally work (with not yet reasonable results). Thole screening to avoid polarization catastrophe. Ions are particularly problematic in this regard. I've seen this mentioned in the Manual, but hadn't ever hit GROMACS specific details on how to apply polarization in the input files. The only source I could locate is within the GROMACS package, under the name sw.itp. It does exclusively implement the so called water polarization model -at least I think so. The water polarization function is a water-specific anisotropy function. Don't try to use it for anything else; the interpretation of the atom numbers for local axis construction are very specific. Can you please direct me to a source from which I could learn how to polarize GROMACS? I was'not lucky on the Internet and, indeed, even at the GROMACS site or Manual (neighter application examples nor file format descriptions). The Thole screening function is (in the released version) not used by anything, so it's not documented. In its present incarnation, you need a [thole_polarization] directive that lists atom-shell/Drude pairs as follows: atom_i shell_i atom_j shell_j 2 a alpha_i alpha_j The 2 is a required function type. My implementation of the CHARMM Drude FF is nearly done, and there are changes to the way the Thole directive is laid out in the future, but at the moment (up through version 5.0), this is the way it works. The code is in src/gromacs/gmxlib/bondfree.c. -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 jalem...@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-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at
Re: [gmx-users] hints for core/shell optimization?
On 7/29/14, 9:33 AM, Tamas Karpati wrote: Dear Justin, I've tried applying the [thole_polarization] section in the format you've suggested. I've noticed the followings. 1. First I made a mistake and forgot to remove the [polarization] section from the TOP file. Even this way MDRUN worked and indicated the usual ca. 10 kJ/mol contribution from polarization. Besides a new contribution of ca. 600 kJ/mol appeared from Thole pol. I switched the former off to have just one kind of pol effect. Is this correct? No, they're not mutually exclusive. The Thole function screens nearby dipoles so they don't experience artificial forces. The thole_polarization name is a bit unfortunate; it should really be thole_screening instead. I suspect you're seeing large forces due to one (or more effects). If your model doesn't have any bonds, then you don't have any excluded interactions, so instead of getting a screened Coulombic interaction via Thole screening, you're getting both added together. You should make use of [exclusions] for neighboring pairs, likely, but I guess the bigger question is how your force field was parametrized and what its expected behavior should be. The other effect is one that I noticed this weekend. I think the forces coming out of the Thole function are incorrect; at least they are in the case of our Drude polarizable force field. I changed the force calculation and everything works as expected now, in my development version of the code. 2. I assumed that in the [thole_polarization] section of the TOP file I needed to select which dipole pairs are to be included on basis of my chemical intuition -if I had any. I tried a few radii of sphere-of-inclusion but saw little differences. Is it a correct selection for Thole screening? (I have a crystal with separate ions held together by Buckingham forces and Thole polarization -there are no bonds at all.) Do I have a way to switch on/off intercell dipole interactions? I don't know what the right answer is here. What is the source of the force field parameters? In our Drude model, 1-2 and 1-3 interactions (based on atom connectivity, not Drude connectivity) are excluded and electrostatics calculated via Thole. You don't have any bonds, so the choice here is not exactly clear and I don't want to pre-judge without knowing what the force field is or how it should be expected to behave. 3. Playing with the a factor of the Thole scheme resulted in more pronounced changes in the contribution from polarization but I still As it should; the value of a should be the sum of the two atomic Thole factors, or 2.6, whichever your force field uses. see too strong Coulombic interactions. Omitting shell particles and thus polarization do scale Coulombics down by a factor of 100 or 1000. I guess that using [polarization] all dipoles interact with all the others (although no idea how dipoles interact through walls of the cell). I imagine that in the [thole_polarization] case only those dipole pairs will interact that I have mindlessly selected. In contrast with my expectation I noticed an order of magnitude greater Thole pol. values than Polarization values in the simple pol. case. (For your information, I have ca. 2000 atoms + 1000 shell particles in a box of ca. 3x3x3 nm. In case of the Thole scheme I have ca. 1e5 pcs. of dipole pairs with a 1 nm radius for inclusion while ca. 15000 for a radius of 0.5 nm.) Sounds like it could be incorrect forces, but without seeing the Thole section of your topology, I can't say for sure. 0.5 - 1.0 nm sounds like far too large of a radius for Thole interactions. -Justin The total potential is almost a pure Coulombic one and is too large. Do you have any hint on what could go wrong? Thanks for your attention. With best regards, toma On Fri, Jul 25, 2014 at 9:31 PM, Tamas Karpati tkarp...@gmail.com wrote: Dear Justin, Thank you for the valuable information. I'm starting experimenting with it. Have a nice weekend, toma On Fri, Jul 25, 2014 at 7:53 PM, Justin Lemkul jalem...@vt.edu wrote: On 7/25/14, 12:40 PM, Tamas Karpati wrote: Dear Justin, Thanks for your educational answers. Coulombic interactions fail at short distances; you probably need to apply I was afraid of that... somehow removing shells from the cores in the initial structure have let it functionally work (with not yet reasonable results). Thole screening to avoid polarization catastrophe. Ions are particularly problematic in this regard. I've seen this mentioned in the Manual, but hadn't ever hit GROMACS specific details on how to apply polarization in the input files. The only source I could locate is within the GROMACS package, under the name sw.itp. It does exclusively implement the so called water polarization model -at least I think so. The water polarization function is a water-specific anisotropy function. Don't try to use it for anything else; the interpretation of
Re: [gmx-users] hints for core/shell optimization?
On 7/25/14, 7:43 AM, Tamas Karpati wrote: Dear all, I have two questions about geometry optimization of a crystal with polarization via the core/shell model. I'm creating *.gro and *.top files by hand and compile them with *.mdp to *.tpr via GROMPP. My FF is also made by hand (simply because i need to learn GROMACS). I have learnt on this list that with Buckingham potentials I need to use the group rather than the Verlet scheme. (1/2) Letting some of all atoms be polarizable through applying shell particles made MDRUN segfault like this: # ... # Reading file AAA_opt.tpr, VERSION 4.6.3 (single precision) # Using 2 MPI threads # # Steepest Descents: # Tolerance (Fmax) = 1.0e+01 #Number of steps= 10 # Segmentation fault1.0e-02 nm, Epot= -nan Fmax= 3.76506e+03, atom= 1357 I imagined some divison by zero situation not handled and have put some random noise on the shell particles' position so they do not anymore start exactly at the atomic sites (meaning nonzero distances). Seemed to work, at least no further crashes. Only energies and forces seem very high: # Steepest Descents: # Tolerance (Fmax) = 1.0e+01 # Number of steps= 10 # Step=0, Dmax= 1.0e-02 nm, Epot= -1.36425e+07 Fmax= 2.99600e+05, atom= 160 # Step=1, Dmax= 1.0e-02 nm, Epot= -1.62080e+07 Fmax= 1.25769e+06, atom= 160 # Step=2, Dmax= 1.2e-02 nm, Epot= -1.95965e+07 Fmax= 6.87820e+08, atom= 2759 # Step=3, Dmax= 1.4e-02 nm, Epot= -2.02902e+07 Fmax= 1.30719e+09, atom= 468 # Step=8, Dmax= 1.1e-03 nm, Epot= -2.18970e+07 Fmax= 5.77722e+09, atom= 1095 # Step= 10, Dmax= 6.5e-04 nm, Epot= -1.96952e+07 Fmax= 3.92889e+08, atom= 1096 # Energy minimization reached the maximum numberof steps before the forces # reached the requestedprecision Fmax 10. My question is the following. Is it (randomized shell positions) a correct procedure with GROMACS? (2/2) Changing from a randomized x/y/z set to a fixed distance at a random direction for the shell particles led to another unexpected result. I scanned a range between 1e-4 to 0.1 nm and noticed that the final core-to-shell distance is a function of the starting one. I used niter = 1 (note: the default is 20) as i noticed in an MD type of job that 20, 100 or 1000 steps were insufficient for the shells to relax within default tolerance. The cell size was ca. 3x3x3 nm. Please provide a full .mdp file; other settings are very relevant. My question is the following. What would be the appropriate core-to-shell distance to apply? The equilibrium distance for core-shell bonds should be zero, deviations from this non-polarized state account for the polarization energy. I appreciate any help so thanks in advance. Upgrade to 4.6.6; there have been issues with shells that have been fixed since 4.6.3. -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 jalem...@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-requ...@gromacs.org.
Re: [gmx-users] hints for core/shell optimization?
Dear Justin, Thank you for your quick answer. Same with GROMACS-4.6.6: core-to-shell distance must be 0 to not crash. My crystal is expected to be polarized (metallic and oxygen sites are the victims of this model). The *.mpd file being used is: # nstcalcenergy = 1 ; 4.6.6 claims this necessary, 4.6.3 didn't need it # integrator = steep # vdw-type= cut-off # coulombtype = pme # nsteps = 10 # periodic-molecules = no # cutoff-scheme = group # ns-type = grid # emtol = 10 ; default=10 kJ/mol/nm # niter = 1 ; default=20 # ;fcstep = ; default=0 ps^2 ; not quite clear what it is # rlist = 1.0 # rcoulomb= 1.0 # rvdw= 1.0 # pbc = xyz # nstxout = 1 ; wouldn't emit pos for each, though # nstfout = 1 ; -- # nstlist = 1000 ; avoid it, no bonds at all # nstlog = 1 # ;pcoupl = no ; used to switch cell-optimization on/off Best regards, toma On Fri, Jul 25, 2014 at 1:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 7/25/14, 7:43 AM, Tamas Karpati wrote: Dear all, I have two questions about geometry optimization of a crystal with polarization via the core/shell model. I'm creating *.gro and *.top files by hand and compile them with *.mdp to *.tpr via GROMPP. My FF is also made by hand (simply because i need to learn GROMACS). I have learnt on this list that with Buckingham potentials I need to use the group rather than the Verlet scheme. (1/2) Letting some of all atoms be polarizable through applying shell particles made MDRUN segfault like this: # ... # Reading file AAA_opt.tpr, VERSION 4.6.3 (single precision) # Using 2 MPI threads # # Steepest Descents: # Tolerance (Fmax) = 1.0e+01 #Number of steps= 10 # Segmentation fault1.0e-02 nm, Epot= -nan Fmax= 3.76506e+03, atom= 1357 I imagined some divison by zero situation not handled and have put some random noise on the shell particles' position so they do not anymore start exactly at the atomic sites (meaning nonzero distances). Seemed to work, at least no further crashes. Only energies and forces seem very high: # Steepest Descents: # Tolerance (Fmax) = 1.0e+01 # Number of steps= 10 # Step=0, Dmax= 1.0e-02 nm, Epot= -1.36425e+07 Fmax= 2.99600e+05, atom= 160 # Step=1, Dmax= 1.0e-02 nm, Epot= -1.62080e+07 Fmax= 1.25769e+06, atom= 160 # Step=2, Dmax= 1.2e-02 nm, Epot= -1.95965e+07 Fmax= 6.87820e+08, atom= 2759 # Step=3, Dmax= 1.4e-02 nm, Epot= -2.02902e+07 Fmax= 1.30719e+09, atom= 468 # Step=8, Dmax= 1.1e-03 nm, Epot= -2.18970e+07 Fmax= 5.77722e+09, atom= 1095 # Step= 10, Dmax= 6.5e-04 nm, Epot= -1.96952e+07 Fmax= 3.92889e+08, atom= 1096 # Energy minimization reached the maximum numberof steps before the forces # reached the requestedprecision Fmax 10. My question is the following. Is it (randomized shell positions) a correct procedure with GROMACS? (2/2) Changing from a randomized x/y/z set to a fixed distance at a random direction for the shell particles led to another unexpected result. I scanned a range between 1e-4 to 0.1 nm and noticed that the final core-to-shell distance is a function of the starting one. I used niter = 1 (note: the default is 20) as i noticed in an MD type of job that 20, 100 or 1000 steps were insufficient for the shells to relax within default tolerance. The cell size was ca. 3x3x3 nm. Please provide a full .mdp file; other settings are very relevant. My question is the following. What would be the appropriate core-to-shell distance to apply? The equilibrium distance for core-shell bonds should be zero, deviations from this non-polarized state account for the polarization energy. I appreciate any help so thanks in advance. Upgrade to 4.6.6; there have been issues with shells that have been fixed since 4.6.3. -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 jalem...@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-requ...@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
Re: [gmx-users] hints for core/shell optimization?
On 7/25/14, 8:10 AM, Tamas Karpati wrote: Dear Justin, Thank you for your quick answer. Same with GROMACS-4.6.6: core-to-shell distance must be 0 to not crash. My crystal is expected to be polarized Does your topology specify the proper intramolecular exclusions? What is(are) the molecule(s)? (metallic and oxygen sites are the victims of this model). The *.mpd file being used is: # nstcalcenergy = 1 ; 4.6.6 claims this necessary, 4.6.3 didn't need it Definitely true. # integrator = steep # vdw-type= cut-off # coulombtype = pme # nsteps = 10 # periodic-molecules = no # cutoff-scheme = group # ns-type = grid # emtol = 10 ; default=10 kJ/mol/nm # niter = 1 ; default=20 # ;fcstep = ; default=0 ps^2 ; not quite clear what it is # rlist = 1.0 # rcoulomb= 1.0 # rvdw= 1.0 # pbc = xyz # nstxout = 1 ; wouldn't emit pos for each, though # nstfout = 1 ; -- # nstlist = 1000 ; avoid it, no bonds at all Try nstlist = 1. The shell positions are solved via SCF (EM), so you need to update the neighbor list very frequently. -Justin # nstlog = 1 # ;pcoupl = no ; used to switch cell-optimization on/off Best regards, toma On Fri, Jul 25, 2014 at 1:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 7/25/14, 7:43 AM, Tamas Karpati wrote: Dear all, I have two questions about geometry optimization of a crystal with polarization via the core/shell model. I'm creating *.gro and *.top files by hand and compile them with *.mdp to *.tpr via GROMPP. My FF is also made by hand (simply because i need to learn GROMACS). I have learnt on this list that with Buckingham potentials I need to use the group rather than the Verlet scheme. (1/2) Letting some of all atoms be polarizable through applying shell particles made MDRUN segfault like this: # ... # Reading file AAA_opt.tpr, VERSION 4.6.3 (single precision) # Using 2 MPI threads # # Steepest Descents: # Tolerance (Fmax) = 1.0e+01 #Number of steps= 10 # Segmentation fault1.0e-02 nm, Epot= -nan Fmax= 3.76506e+03, atom= 1357 I imagined some divison by zero situation not handled and have put some random noise on the shell particles' position so they do not anymore start exactly at the atomic sites (meaning nonzero distances). Seemed to work, at least no further crashes. Only energies and forces seem very high: # Steepest Descents: # Tolerance (Fmax) = 1.0e+01 # Number of steps= 10 # Step=0, Dmax= 1.0e-02 nm, Epot= -1.36425e+07 Fmax= 2.99600e+05, atom= 160 # Step=1, Dmax= 1.0e-02 nm, Epot= -1.62080e+07 Fmax= 1.25769e+06, atom= 160 # Step=2, Dmax= 1.2e-02 nm, Epot= -1.95965e+07 Fmax= 6.87820e+08, atom= 2759 # Step=3, Dmax= 1.4e-02 nm, Epot= -2.02902e+07 Fmax= 1.30719e+09, atom= 468 # Step=8, Dmax= 1.1e-03 nm, Epot= -2.18970e+07 Fmax= 5.77722e+09, atom= 1095 # Step= 10, Dmax= 6.5e-04 nm, Epot= -1.96952e+07 Fmax= 3.92889e+08, atom= 1096 # Energy minimization reached the maximum numberof steps before the forces # reached the requestedprecision Fmax 10. My question is the following. Is it (randomized shell positions) a correct procedure with GROMACS? (2/2) Changing from a randomized x/y/z set to a fixed distance at a random direction for the shell particles led to another unexpected result. I scanned a range between 1e-4 to 0.1 nm and noticed that the final core-to-shell distance is a function of the starting one. I used niter = 1 (note: the default is 20) as i noticed in an MD type of job that 20, 100 or 1000 steps were insufficient for the shells to relax within default tolerance. The cell size was ca. 3x3x3 nm. Please provide a full .mdp file; other settings are very relevant. My question is the following. What would be the appropriate core-to-shell distance to apply? The equilibrium distance for core-shell bonds should be zero, deviations from this non-polarized state account for the polarization energy. I appreciate any help so thanks in advance. Upgrade to 4.6.6; there have been issues with shells that have been fixed since 4.6.3. -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 jalem...@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
Re: [gmx-users] hints for core/shell optimization?
On 7/25/14, 11:22 AM, Tamas Karpati wrote: Dear Justin, Does your topology specify the proper intramolecular exclusions? What is(are) the molecule(s)? No bonds, no exclusions. The whole crystal is modelled by ions interacting via forces of the Coulomb and Buckingham types. In fact, there is an X-Y-X angle force type which does have an effect on forces and energy. Shell particles, on the other hand, are defined via the [polarization] section as quasi-bonds. That's all. Coulombic interactions fail at short distances; you probably need to apply Thole screening to avoid polarization catastrophe. Ions are particularly problematic in this regard. Try nstlist = 1. The shell positions are solved via SCF (EM), so you need to update the neighbor list very frequently. Thanks for the trick. Tried but to no avail (very same results). Although shells should have their weight in the model, I expect them not to frequently change partner -each is ultimately connected to its core atomic site. In addition, there are no more bonds in the system i, thus, don't see the point in regenerating pairs. Is there a process -I should get to know about- which autogenerates bonds? If your model doesn't use them, then there's nothing to be done here. I was asking about bonded structure and exclusions and such because of the Thole issue I noted above. -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 jalem...@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-requ...@gromacs.org.
Re: [gmx-users] hints for core/shell optimization?
Dear Justin, Does your topology specify the proper intramolecular exclusions? What is(are) the molecule(s)? No bonds, no exclusions. The whole crystal is modelled by ions interacting via forces of the Coulomb and Buckingham types. In fact, there is an X-Y-X angle force type which does have an effect on forces and energy. Shell particles, on the other hand, are defined via the [polarization] section as quasi-bonds. That's all. Try nstlist = 1. The shell positions are solved via SCF (EM), so you need to update the neighbor list very frequently. Thanks for the trick. Tried but to no avail (very same results). Although shells should have their weight in the model, I expect them not to frequently change partner -each is ultimately connected to its core atomic site. In addition, there are no more bonds in the system i, thus, don't see the point in regenerating pairs. Is there a process -I should get to know about- which autogenerates bonds? Thanks al lot for you help. Best regards, toma -- 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-requ...@gromacs.org.
Re: [gmx-users] hints for core/shell optimization?
Dear Justin, Thanks for your educational answers. Coulombic interactions fail at short distances; you probably need to apply I was afraid of that... somehow removing shells from the cores in the initial structure have let it functionally work (with not yet reasonable results). Thole screening to avoid polarization catastrophe. Ions are particularly problematic in this regard. I've seen this mentioned in the Manual, but hadn't ever hit GROMACS specific details on how to apply polarization in the input files. The only source I could locate is within the GROMACS package, under the name sw.itp. It does exclusively implement the so called water polarization model -at least I think so. Can you please direct me to a source from which I could learn how to polarize GROMACS? I was'not lucky on the Internet and, indeed, even at the GROMACS site or Manual (neighter application examples nor file format descriptions). I appreciate your help so much. With regards, toma -- 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-requ...@gromacs.org.
Re: [gmx-users] hints for core/shell optimization?
On 7/25/14, 12:40 PM, Tamas Karpati wrote: Dear Justin, Thanks for your educational answers. Coulombic interactions fail at short distances; you probably need to apply I was afraid of that... somehow removing shells from the cores in the initial structure have let it functionally work (with not yet reasonable results). Thole screening to avoid polarization catastrophe. Ions are particularly problematic in this regard. I've seen this mentioned in the Manual, but hadn't ever hit GROMACS specific details on how to apply polarization in the input files. The only source I could locate is within the GROMACS package, under the name sw.itp. It does exclusively implement the so called water polarization model -at least I think so. The water polarization function is a water-specific anisotropy function. Don't try to use it for anything else; the interpretation of the atom numbers for local axis construction are very specific. Can you please direct me to a source from which I could learn how to polarize GROMACS? I was'not lucky on the Internet and, indeed, even at the GROMACS site or Manual (neighter application examples nor file format descriptions). The Thole screening function is (in the released version) not used by anything, so it's not documented. In its present incarnation, you need a [thole_polarization] directive that lists atom-shell/Drude pairs as follows: atom_i shell_i atom_j shell_j 2 a alpha_i alpha_j The 2 is a required function type. My implementation of the CHARMM Drude FF is nearly done, and there are changes to the way the Thole directive is laid out in the future, but at the moment (up through version 5.0), this is the way it works. The code is in src/gromacs/gmxlib/bondfree.c. -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 jalem...@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-requ...@gromacs.org.
