From a4a97de84f561c9b65689bd477298198bd0f4e81 Mon Sep 17 00:00:00 2001
From: Michele Ceriotti <michele.ceriotti@gmail.com>
Date: Wed, 8 Apr 2015 10:45:49 +0200
Subject: [PATCH] A few GLE fixes
---
doc/fix_gle.html | 16 +++++++++++++++-
doc/fix_gle.txt | 16 +++++++++++++++-
src/USER-MISC/fix_gle.cpp | 8 +++++++-
src/USER-MISC/fix_ipi.cpp | 14 ++++++--------
4 files changed, 43 insertions(+), 11 deletions(-)
diff --git a/doc/fix_gle.html b/doc/fix_gle.html
index e82902a973..5ae0c8783d 100644
--- a/doc/fix_gle.html
+++ b/doc/fix_gle.html
@@ -50,7 +50,8 @@ vibrational modes in the system to inexpensive (approximate)
modelling of nuclear quantum effects. Contrary to
<A HREF = "fix_langevin.html">fix langevin</A>, this fix performs both
thermostatting and evolution of the Hamiltonian equations of motion, so it
-does not need to be used together with <A HREF = "fix_nve.html">fix nve</A>.
+should not be used together with <A HREF = "fix_nve.html">fix nve</A> -- at least not
+on the same atom groups.
</P>
<P>Each degree of freedom in the thermostatted group is supplemented
with Ns additional degrees of freedom s, and the equations of motion
@@ -93,6 +94,14 @@ generate its own unique seed and its own stream of random numbers.
Thus the dynamics of the system will not be identical on two runs on
different numbers of processors.
</P>
+<P>Note also that the Generalized Langevin Dynamics scheme that is
+implemented by the <A HREF = "fix_gld.html">fix gld</A> scheme is closely related
+to the present one. In fact, it should be always possible to cast the
+Prony series form of the memory kernel used by GLD into an appropriate
+input matrix for <A HREF = "fix_gle.html">fix_gle</A>. While the GLE scheme is more
+general, the form used by <A HREF = "fix_gld.html">fix gld</A> can be more directly
+related to the representation of an implicit solvent environment.
+</P>
<P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
</P>
<P>The instantaneous values of the extended variables are written to
@@ -127,6 +136,11 @@ with rigid bodies, SHAKE or RATTLE. It is expected that all the
thermostatted degrees of freedom are fully flexible, and the sampled
ensemble will not be correct otherwise.
</P>
+<P>In order to perform constant-pressure simulations please use
+<A HREF = "fix_press_berendsen.html">fix press/berendsen</A>, rather than
+<A HREF = "fix_npt.html">fix_npt</A>, to avoid duplicate integration of the
+equations of motion.
+</P>
<P>This fix is part of the USER-MISC package. It is only enabled if LAMMPS
was built with that package. See the <A HREF = "Section_start.html#start_3">Making
LAMMPS</A> section for more info.
diff --git a/doc/fix_gle.txt b/doc/fix_gle.txt
index 08ff64777c..8ce2a2b6fd 100644
--- a/doc/fix_gle.txt
+++ b/doc/fix_gle.txt
@@ -40,7 +40,8 @@ vibrational modes in the system to inexpensive (approximate)
modelling of nuclear quantum effects. Contrary to
"fix langevin"_fix_langevin.html, this fix performs both
thermostatting and evolution of the Hamiltonian equations of motion, so it
-does not need to be used together with "fix nve"_fix_nve.html.
+should not be used together with "fix nve"_fix_nve.html -- at least not
+on the same atom groups.
Each degree of freedom in the thermostatted group is supplemented
with Ns additional degrees of freedom s, and the equations of motion
@@ -83,6 +84,14 @@ generate its own unique seed and its own stream of random numbers.
Thus the dynamics of the system will not be identical on two runs on
different numbers of processors.
+Note also that the Generalized Langevin Dynamics scheme that is
+implemented by the "fix gld"_fix_gld.html scheme is closely related
+to the present one. In fact, it should be always possible to cast the
+Prony series form of the memory kernel used by GLD into an appropriate
+input matrix for "fix_gle"_fix_gle.html. While the GLE scheme is more
+general, the form used by "fix gld"_fix_gld.html can be more directly
+related to the representation of an implicit solvent environment.
+
[Restart, fix_modify, output, run start/stop, minimize info:]
The instantaneous values of the extended variables are written to
@@ -117,6 +126,11 @@ with rigid bodies, SHAKE or RATTLE. It is expected that all the
thermostatted degrees of freedom are fully flexible, and the sampled
ensemble will not be correct otherwise.
+In order to perform constant-pressure simulations please use
+"fix press/berendsen"_fix_press_berendsen.html, rather than
+"fix_npt"_fix_npt.html, to avoid duplicate integration of the
+equations of motion.
+
This fix is part of the USER-MISC package. It is only enabled if LAMMPS
was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
diff --git a/src/USER-MISC/fix_gle.cpp b/src/USER-MISC/fix_gle.cpp
index 5f28dea37b..95f868f95f 100644
--- a/src/USER-MISC/fix_gle.cpp
+++ b/src/USER-MISC/fix_gle.cpp
@@ -82,7 +82,13 @@ void StabCholesky(int n, const double* MMt, double* M)
for (k=0; k<i; ++k) D[i]-=L[midx(n,i,k)]*L[midx(n,i,k)]*D[k];
}
- for(i=0; i<n; ++i) D[i]=(D[i]>0.?sqrt(D[i]):0.);
+ for(i=0; i<n; ++i)
+ {
+#ifdef GLE_DEBUG
+ if (D[i]<0) fprintf(stderr,"GLE Cholesky: Negative diagonal term %le, has been set to zero.\n", D[i]);
+#endif
+ D[i]=(D[i]>0.?sqrt(D[i]):0.);
+ }
for(i=0; i<n; ++i) for (j=0; j<n; j++) M[midx(n,i,j)]=L[midx(n,i,j)]*D[j];
}
diff --git a/src/USER-MISC/fix_ipi.cpp b/src/USER-MISC/fix_ipi.cpp
index 02cfcd9abb..58251a814e 100644
--- a/src/USER-MISC/fix_ipi.cpp
+++ b/src/USER-MISC/fix_ipi.cpp
@@ -319,12 +319,12 @@ void FixIPI::initial_integrate(int vflag)
double *boxlo = domain->boxlo;
double posconv;
posconv=0.52917721*force->angstrom;
- boxlo[0] = 0;
- boxlo[1] = 0;
- boxlo[2] = 0;
- boxhi[0] = cellh[0]*posconv;
- boxhi[1] = cellh[4]*posconv;
- boxhi[2] = cellh[8]*posconv;
+ boxlo[0] = -0.5*cellh[0]*posconv;
+ boxlo[1] = -0.5*cellh[4]*posconv;
+ boxlo[2] = -0.5*cellh[8]*posconv;
+ boxhi[0] = -boxlo[0];
+ boxhi[1] = -boxlo[1];
+ boxhi[2] = -boxlo[2];
domain->xy = cellh[1]*posconv;
domain->xz = cellh[2]*posconv;
domain->yz = cellh[5]*posconv;
@@ -421,7 +421,6 @@ void FixIPI::final_integrate()
error->one(FLERR, "Got EXIT message from i-PI. Now leaving!");
if (strcmp(header,"GETFORCE ") == 0 ) {
-
writebuffer(ipisock,"FORCEREADY ",MSGLEN, error);
writebuffer(ipisock,(char*) &pot,8, error);
writebuffer(ipisock,(char*) &nat,4, error);
@@ -432,7 +431,6 @@ void FixIPI::final_integrate()
}
else
error->one(FLERR, "Wrapper did not ask for forces, I will now die!");
-
}
hasdata=0;
--
GitLab