git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12775 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
@ -36,7 +36,9 @@ dihedral_style quadratic, Loukas Peristeras, loukas.peristeras at scienomics.com
|
|||||||
dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12
|
dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12
|
||||||
fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
|
fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
|
||||||
fix ave/spatial/sphere, Niall Jackson (Imperial College), niall.jackson at gmail.com, 18 Nov 14
|
fix ave/spatial/sphere, Niall Jackson (Imperial College), niall.jackson at gmail.com, 18 Nov 14
|
||||||
|
fix gle, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014
|
||||||
fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009
|
fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009
|
||||||
|
fix ipi, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014
|
||||||
fix smd, Axel Kohlmeyer, akohlmey at gmail.com, 19 May 2008
|
fix smd, Axel Kohlmeyer, akohlmey at gmail.com, 19 May 2008
|
||||||
fix ti/rs, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013
|
fix ti/rs, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013
|
||||||
fix ti/spring, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013
|
fix ti/spring, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013
|
||||||
|
|||||||
868
src/USER-MISC/fix_gle.cpp
Normal file
868
src/USER-MISC/fix_gle.cpp
Normal file
@ -0,0 +1,868 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||||
|
http://lammps.sandia.gov, Sandia National Laboratories
|
||||||
|
Steve Plimpton, sjplimp@sandia.gov
|
||||||
|
|
||||||
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||||
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||||
|
certain rights in this software. This software is distributed under
|
||||||
|
the GNU General Public License.
|
||||||
|
|
||||||
|
See the README file in the top-level LAMMPS directory.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Contributing authors: Michele Ceriotti (EPFL), Joe Morrone (Stony Brook)
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#include "mpi.h"
|
||||||
|
#include "math.h"
|
||||||
|
#include "string.h"
|
||||||
|
#include "stdlib.h"
|
||||||
|
#include "fix_gle.h"
|
||||||
|
#include "math_extra.h"
|
||||||
|
#include "atom.h"
|
||||||
|
#include "atom_vec_ellipsoid.h"
|
||||||
|
#include "force.h"
|
||||||
|
#include "update.h"
|
||||||
|
#include "modify.h"
|
||||||
|
#include "compute.h"
|
||||||
|
#include "domain.h"
|
||||||
|
#include "region.h"
|
||||||
|
#include "respa.h"
|
||||||
|
#include "comm.h"
|
||||||
|
#include "input.h"
|
||||||
|
#include "variable.h"
|
||||||
|
#include "random_mars.h"
|
||||||
|
#include "memory.h"
|
||||||
|
#include "error.h"
|
||||||
|
#include "group.h"
|
||||||
|
|
||||||
|
using namespace LAMMPS_NS;
|
||||||
|
using namespace FixConst;
|
||||||
|
|
||||||
|
enum{NOBIAS,BIAS};
|
||||||
|
enum{CONSTANT,EQUAL,ATOM};
|
||||||
|
//#define GLE_DEBUG 1
|
||||||
|
|
||||||
|
#define MAXLINE 1024
|
||||||
|
|
||||||
|
/* syntax for fix_gle:
|
||||||
|
* fix nfix id-group gle ns Tstart Tstop seed amatrix [noneq cmatrix] [every nmts]
|
||||||
|
*
|
||||||
|
* */
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
namespace GLE {
|
||||||
|
/* GLE Utility functions -- might perhaps use standard LAMMPS routines if available */
|
||||||
|
#define midx(n,i,j) ((i)*(n)+(j))
|
||||||
|
|
||||||
|
//"stabilized" cholesky decomposition. does a LDL^t decomposition, then sets to zero the negative diagonal elements and gets MM^t
|
||||||
|
void StabCholesky(int n, const double* MMt, double* M)
|
||||||
|
{
|
||||||
|
double L[n*n], D[n];
|
||||||
|
|
||||||
|
int i,j,k;
|
||||||
|
for(i=0; i<n; ++i) D[i]=0.;
|
||||||
|
for(i=0; i<n*n; ++i) L[i]=0.;
|
||||||
|
|
||||||
|
for(i=0; i<n; ++i)
|
||||||
|
{
|
||||||
|
L[midx(n,i,i)]=1.;
|
||||||
|
for (j=0; j<i; j++)
|
||||||
|
{
|
||||||
|
L[midx(n,i,j)]=MMt[midx(n,i,j)];
|
||||||
|
for (k=0; k<j; ++k) L[midx(n,i,j)]-=L[midx(n,i,k)]*L[midx(n,j,k)]*D[k];
|
||||||
|
if (D[j]!=0.) L[midx(n,i,j)]/=D[j];
|
||||||
|
else L[midx(n,i,j)]=0.0;
|
||||||
|
}
|
||||||
|
D[i]=MMt[midx(n,i,i)];
|
||||||
|
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) for (j=0; j<n; j++) M[midx(n,i,j)]=L[midx(n,i,j)]*D[j];
|
||||||
|
}
|
||||||
|
|
||||||
|
void MyMult(int n, int m, int r, const double* A, const double* B, double* C, double cf=0.0)
|
||||||
|
{
|
||||||
|
// !! TODO !! should probably call BLAS (or check if some efficient matrix-matrix multiply is implemented somewhere in LAMMPS)
|
||||||
|
// (rows x cols) :: A is n x r, B is r x m, C is n x m
|
||||||
|
int i,j,k; double *cij;
|
||||||
|
for (i=0; i<n; ++i)
|
||||||
|
for (j=0; j<m; ++j)
|
||||||
|
{
|
||||||
|
cij=&C[midx(m,i,j)]; *cij *= cf;
|
||||||
|
for (k=0; k<r; ++k) *cij+=A[midx(r,i,k)]*B[midx(m,k,j)];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
#define MIN(A,B) ((A) < (B) ? (A) : (B))
|
||||||
|
// BLAS-like version of MyMult(). AK 2014-08-06
|
||||||
|
|
||||||
|
inline void AkMult(const int n, const int m, const int r,
|
||||||
|
const double * const A, const double * const B,
|
||||||
|
double * const C, const double cf=0.0)
|
||||||
|
{
|
||||||
|
// block buffer
|
||||||
|
const int blk=64;
|
||||||
|
double buf[blk*blk];
|
||||||
|
int i,j,k,ib,jb,kb;
|
||||||
|
|
||||||
|
for (i = 0; i < (n*m); ++i) C[i] *= cf;
|
||||||
|
|
||||||
|
for (kb = 0; kb < r; kb +=blk) {
|
||||||
|
// fill buffer
|
||||||
|
const int ke = MIN(kb+blk,r);
|
||||||
|
for (ib = 0; ib < n; ib += blk) {
|
||||||
|
const int ie = MIN(ib+blk,n);
|
||||||
|
for (i = ib; i < ie; ++i) {
|
||||||
|
for (k = kb; k < ke; ++k) {
|
||||||
|
buf[midx(blk,k-kb,i-ib)] = A[midx(r,i,k)];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
for (jb = 0; jb < m; jb += blk) {
|
||||||
|
double tmp;
|
||||||
|
const int je = MIN(jb+blk,m);
|
||||||
|
|
||||||
|
for (j = jb; j < je; ++j) {
|
||||||
|
for (i = ib; i < ie; ++i) {
|
||||||
|
tmp = 0.0;
|
||||||
|
for (k = kb; k < ke; ++k)
|
||||||
|
tmp += buf[midx(blk,k-kb,i-ib)] * B[midx(m,k,j)];
|
||||||
|
|
||||||
|
C[midx(m,i,j)] += tmp;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void MyTrans(int n, const double* A, double* AT)
|
||||||
|
{
|
||||||
|
for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) AT[j*n+i]=A[i*n+j];
|
||||||
|
}
|
||||||
|
|
||||||
|
void MyPrint(int n, const double* A)
|
||||||
|
{
|
||||||
|
for (int k=0; k<n*n; ++k) { printf("%10.5e ", A[k]); if ((k+1)%n==0) printf("\n");}
|
||||||
|
}
|
||||||
|
|
||||||
|
//matrix exponential by scaling and squaring.
|
||||||
|
void MatrixExp(int n, const double* M, double* EM, int j=8, int k=8)
|
||||||
|
{
|
||||||
|
double tc[j+1], SM[n*n], TMP[n*n];
|
||||||
|
double onetotwok=pow(0.5,1.0*k);
|
||||||
|
|
||||||
|
|
||||||
|
tc[0]=1;
|
||||||
|
for (int i=0; i<j; ++i) tc[i+1]=tc[i]/(i+1.0);
|
||||||
|
|
||||||
|
for (int i=0; i<n*n; ++i) { SM[i]=M[i]*onetotwok; EM[i]=0.0; TMP[i]=0.0; }
|
||||||
|
|
||||||
|
for (int i=0; i<n; ++i) EM[midx(n,i,i)]=tc[j];
|
||||||
|
|
||||||
|
//taylor exp of scaled matrix
|
||||||
|
for (int p=j-1; p>=0; p--)
|
||||||
|
{
|
||||||
|
MyMult(n, n, n, SM, EM, TMP); for (int i=0; i<n*n; ++i) EM[i]=TMP[i];
|
||||||
|
for (int i=0; i<n; ++i) EM[midx(n,i,i)]+=tc[p];
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int p=0; p<k; p++)
|
||||||
|
{ MyMult(n, n, n, EM, EM, TMP); for (int i=0; i<n*n; ++i) EM[i]=TMP[i]; }
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) :
|
||||||
|
Fix(lmp, narg, arg)
|
||||||
|
{
|
||||||
|
if (narg < 8)
|
||||||
|
error->all(FLERR,"Illegal fix gle command. Expecting: fix <fix-ID>"
|
||||||
|
" <group-ID> gle <ns> <Tstart> <Tstop> <seed> <Amatrix>");
|
||||||
|
|
||||||
|
restart_peratom = 1;
|
||||||
|
time_integrate = 1;
|
||||||
|
|
||||||
|
// number of additional momenta
|
||||||
|
ns = force->inumeric(FLERR,arg[3]);
|
||||||
|
ns1sq = (ns+1)*(ns+1);
|
||||||
|
|
||||||
|
// allocate GLE matrices
|
||||||
|
A = new double[ns1sq];
|
||||||
|
C = new double[ns1sq];
|
||||||
|
T = new double[ns1sq];
|
||||||
|
S = new double[ns1sq];
|
||||||
|
TT = new double[ns1sq];
|
||||||
|
ST = new double[ns1sq];
|
||||||
|
|
||||||
|
// start temperature (t ramp)
|
||||||
|
t_start = force->numeric(FLERR,arg[4]);
|
||||||
|
|
||||||
|
// final temperature (t ramp)
|
||||||
|
t_stop = force->numeric(FLERR,arg[5]);
|
||||||
|
|
||||||
|
// PRNG seed
|
||||||
|
int seed = force->inumeric(FLERR,arg[6]);
|
||||||
|
|
||||||
|
// LOADING A matrix
|
||||||
|
FILE *fgle = NULL;
|
||||||
|
char *fname = arg[7];
|
||||||
|
if (comm->me == 0) {
|
||||||
|
fgle = force->open_potential(fname);
|
||||||
|
if (fgle == NULL) {
|
||||||
|
char str[128];
|
||||||
|
sprintf(str,"Cannot open A-matrix file %s",fname);
|
||||||
|
error->one(FLERR,str);
|
||||||
|
}
|
||||||
|
if (screen) fprintf(screen,"Reading A-matrix from %s\n", fname);
|
||||||
|
if (logfile) fprintf(logfile,"Reading A-matrix from %s\n", fname);
|
||||||
|
}
|
||||||
|
|
||||||
|
// read each line of the file, skipping blank lines or leading '#'
|
||||||
|
|
||||||
|
char line[MAXLINE],*ptr;
|
||||||
|
int n,nwords,ndone=0,eof=0;
|
||||||
|
while (1) {
|
||||||
|
if (comm->me == 0) {
|
||||||
|
ptr = fgets(line,MAXLINE,fgle);
|
||||||
|
if (ptr == NULL) {
|
||||||
|
eof = 1;
|
||||||
|
fclose(fgle);
|
||||||
|
} else n = strlen(line) + 1;
|
||||||
|
}
|
||||||
|
MPI_Bcast(&eof,1,MPI_INT,0,world);
|
||||||
|
if (eof) break;
|
||||||
|
MPI_Bcast(&n,1,MPI_INT,0,world);
|
||||||
|
MPI_Bcast(line,n,MPI_CHAR,0,world);
|
||||||
|
|
||||||
|
// strip comment, skip line if blank
|
||||||
|
|
||||||
|
if ((ptr = strchr(line,'#'))) *ptr = '\0';
|
||||||
|
|
||||||
|
nwords = atom->count_words(line);
|
||||||
|
if (nwords == 0) continue;
|
||||||
|
|
||||||
|
ptr = strtok(line," \t\n\r\f");
|
||||||
|
do {
|
||||||
|
A[ndone] = atof(ptr);
|
||||||
|
ptr = strtok(NULL," \t\n\r\f");
|
||||||
|
ndone++;
|
||||||
|
} while ((ptr != NULL) && (ndone < ns1sq));
|
||||||
|
}
|
||||||
|
|
||||||
|
fnoneq=0; gle_every=1; gle_step=0;
|
||||||
|
for (int iarg=8; iarg<narg; iarg+=2) {
|
||||||
|
if(strcmp(arg[iarg],"noneq") == 0) {
|
||||||
|
fnoneq = 1;
|
||||||
|
if (iarg+2>narg)
|
||||||
|
error->all(FLERR,"Did not specify C matrix for non-equilibrium GLE");
|
||||||
|
|
||||||
|
fname = arg[iarg+1];
|
||||||
|
|
||||||
|
} else if (strcmp(arg[iarg],"every") == 0) {
|
||||||
|
|
||||||
|
if (iarg+2>narg)
|
||||||
|
error->all(FLERR, "Did not specify interval for applying the GLE");
|
||||||
|
gle_every=force->inumeric(FLERR,arg[iarg+1]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// set C matrix
|
||||||
|
|
||||||
|
if (fnoneq == 0) {
|
||||||
|
t_target=t_start;
|
||||||
|
const double kT = t_target * force->boltz / force->mvv2e;
|
||||||
|
memset(C,0,sizeof(double)*ns1sq);
|
||||||
|
for (int i=0; i<ns1sq; i+=(ns+2))
|
||||||
|
C[i]=kT;
|
||||||
|
|
||||||
|
} else {
|
||||||
|
if (comm->me == 0) {
|
||||||
|
fgle = force->open_potential(fname);
|
||||||
|
if (fgle == NULL) {
|
||||||
|
char str[128];
|
||||||
|
sprintf(str,"Cannot open C-matrix file %s",fname);
|
||||||
|
error->one(FLERR,str);
|
||||||
|
}
|
||||||
|
if (screen)
|
||||||
|
fprintf(screen,"Reading C-matrix from %s\n", fname);
|
||||||
|
if (logfile)
|
||||||
|
fprintf(logfile,"Reading C-matrix from %s\n", fname);
|
||||||
|
}
|
||||||
|
|
||||||
|
// read each line of the file, skipping blank lines or leading '#'
|
||||||
|
ndone = eof = 0;
|
||||||
|
const double cfac = force->boltz / force->mvv2e;
|
||||||
|
|
||||||
|
while (1) {
|
||||||
|
if (comm->me == 0) {
|
||||||
|
ptr = fgets(line,MAXLINE,fgle);
|
||||||
|
if (ptr == NULL) {
|
||||||
|
eof = 1;
|
||||||
|
fclose(fgle);
|
||||||
|
} else n = strlen(line) + 1;
|
||||||
|
}
|
||||||
|
MPI_Bcast(&eof,1,MPI_INT,0,world);
|
||||||
|
if (eof) break;
|
||||||
|
MPI_Bcast(&n,1,MPI_INT,0,world);
|
||||||
|
MPI_Bcast(line,n,MPI_CHAR,0,world);
|
||||||
|
|
||||||
|
// strip comment, skip line if blank
|
||||||
|
|
||||||
|
if ((ptr = strchr(line,'#'))) *ptr = '\0';
|
||||||
|
|
||||||
|
nwords = atom->count_words(line);
|
||||||
|
if (nwords == 0) continue;
|
||||||
|
|
||||||
|
ptr = strtok(line," \t\n\r\f");
|
||||||
|
do {
|
||||||
|
C[ndone] = cfac*atof(ptr);
|
||||||
|
ptr = strtok(NULL," \t\n\r\f");
|
||||||
|
ndone++;
|
||||||
|
} while ((ptr != NULL) && (ndone < ns1sq));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
#ifdef GLE_DEBUG
|
||||||
|
printf("A Matrix\n");
|
||||||
|
GLE::MyPrint(ns+1,A);
|
||||||
|
printf("C Matrix\n");
|
||||||
|
GLE::MyPrint(ns+1,C);
|
||||||
|
#endif
|
||||||
|
|
||||||
|
// initialize Marsaglia RNG with processor-unique seed
|
||||||
|
// NB: this means runs will not be the same with different numbers of processors
|
||||||
|
if (seed <= 0) error->all(FLERR,"Illegal fix gle command");
|
||||||
|
random = new RanMars(lmp,seed + comm->me);
|
||||||
|
|
||||||
|
// allocate per-type arrays for mass-scaling
|
||||||
|
sqrt_m=NULL;
|
||||||
|
memory->grow(sqrt_m, atom->ntypes+1,"gle:sqrt_m");
|
||||||
|
|
||||||
|
// allocates space for additional degrees of freedom
|
||||||
|
gle_s=NULL;
|
||||||
|
// allocates space for temporaries
|
||||||
|
gle_tmp1=gle_tmp2=NULL;
|
||||||
|
|
||||||
|
grow_arrays(atom->nmax);
|
||||||
|
init_gles();
|
||||||
|
|
||||||
|
// add callbacks to enable restarts
|
||||||
|
atom->add_callback(0);
|
||||||
|
atom->add_callback(1);
|
||||||
|
|
||||||
|
energy = 0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* --- Frees up memory used by temporaries and buffers ------------------ */
|
||||||
|
|
||||||
|
FixGLE::~FixGLE()
|
||||||
|
{
|
||||||
|
delete random;
|
||||||
|
delete [] A;
|
||||||
|
delete [] C;
|
||||||
|
delete [] S;
|
||||||
|
delete [] T;
|
||||||
|
delete [] TT;
|
||||||
|
delete [] ST;
|
||||||
|
|
||||||
|
memory->destroy(sqrt_m);
|
||||||
|
memory->destroy(gle_s);
|
||||||
|
memory->destroy(gle_tmp1);
|
||||||
|
memory->destroy(gle_tmp2);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixGLE::setmask()
|
||||||
|
{
|
||||||
|
int mask = 0;
|
||||||
|
mask |= INITIAL_INTEGRATE;
|
||||||
|
mask |= FINAL_INTEGRATE;
|
||||||
|
mask |= INITIAL_INTEGRATE_RESPA;
|
||||||
|
mask |= FINAL_INTEGRATE_RESPA;
|
||||||
|
mask |= THERMO_ENERGY;
|
||||||
|
|
||||||
|
|
||||||
|
return mask;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ------- Initializes one-time quantities for GLE ---------------------- */
|
||||||
|
|
||||||
|
void FixGLE::init()
|
||||||
|
{
|
||||||
|
|
||||||
|
dogle = 1;
|
||||||
|
dtv = update->dt;
|
||||||
|
dtf = 0.5 * update->dt * force->ftm2v;
|
||||||
|
|
||||||
|
// set force prefactors
|
||||||
|
if (!atom->rmass) {
|
||||||
|
for (int i = 1; i <= atom->ntypes; i++) {
|
||||||
|
sqrt_m[i] = sqrt(atom->mass[i]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if (strstr(update->integrate_style,"respa")) {
|
||||||
|
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
||||||
|
step_respa = ((Respa *) update->integrate)->step;
|
||||||
|
}
|
||||||
|
|
||||||
|
init_gle();
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ------- Initializes integrator matrices (change with T and dt) ------- */
|
||||||
|
|
||||||
|
void FixGLE::init_gle()
|
||||||
|
{
|
||||||
|
// compute Langevin terms
|
||||||
|
|
||||||
|
double *tmp1 = new double[ns1sq];
|
||||||
|
double *tmp2 = new double[ns1sq];
|
||||||
|
|
||||||
|
for (int i=0; i<ns1sq; ++i) {
|
||||||
|
tmp1[i]=-A[i]*update->dt*0.5*gle_every;
|
||||||
|
tmp2[i]=S[i]=0.0;
|
||||||
|
}
|
||||||
|
GLE::MatrixExp(ns+1,tmp1,T);
|
||||||
|
|
||||||
|
GLE::MyMult(ns+1,ns+1,ns+1,T,C,tmp1);
|
||||||
|
GLE::MyTrans(ns+1,T,tmp2);
|
||||||
|
GLE::MyMult(ns+1,ns+1,ns+1,tmp1,tmp2,S);
|
||||||
|
|
||||||
|
for (int i=0; i<ns1sq; ++i) tmp1[i]=C[i]-S[i];
|
||||||
|
|
||||||
|
GLE::StabCholesky(ns+1, tmp1, S); //!TODO use symmetric square root, which is more stable.
|
||||||
|
|
||||||
|
#ifdef GLE_DEBUG
|
||||||
|
printf("T Matrix\n");
|
||||||
|
GLE::MyPrint(ns+1,T);
|
||||||
|
printf("S Matrix\n");
|
||||||
|
GLE::MyPrint(ns+1,S);
|
||||||
|
#endif
|
||||||
|
|
||||||
|
// transposed evolution matrices to have fast index multiplication in gle_integrate
|
||||||
|
GLE::MyTrans(ns+1,T,TT);
|
||||||
|
GLE::MyTrans(ns+1,S,ST);
|
||||||
|
delete[] tmp1;
|
||||||
|
delete[] tmp2;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ------- Sets initial values of additional DOF for free particles ----- */
|
||||||
|
void FixGLE::init_gles()
|
||||||
|
{
|
||||||
|
|
||||||
|
int *mask = atom->mask;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
double *rootC = new double[ns1sq];
|
||||||
|
double *rootCT = new double[ns1sq];
|
||||||
|
double *newg = new double[3*(ns+1)*nlocal];
|
||||||
|
double *news = new double[3*(ns+1)*nlocal];
|
||||||
|
|
||||||
|
GLE::StabCholesky(ns+1, C, rootC);
|
||||||
|
GLE::MyTrans(ns+1,rootC,rootCT);
|
||||||
|
|
||||||
|
memset(news,0,sizeof(double)*3*(ns+1)*nlocal);
|
||||||
|
for (int i = 0; i < nlocal*3*(ns+1); ++i)
|
||||||
|
newg[i] = random->gaussian();
|
||||||
|
|
||||||
|
GLE::AkMult(nlocal*3,ns+1,ns+1, newg, rootCT, news);
|
||||||
|
|
||||||
|
int nk=0; // unpacks temporary into gle_s
|
||||||
|
for (int i = 0; i < nlocal; i++) {
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
// first loads velocities
|
||||||
|
for (int k = 0; k<3; k++) {
|
||||||
|
for (int j=0; j<ns; ++j)
|
||||||
|
gle_s[i][k*ns+j]=news[nk++];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
delete[] rootC;
|
||||||
|
delete[] rootCT;
|
||||||
|
delete[] news;
|
||||||
|
delete[] newg;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixGLE::setup(int vflag)
|
||||||
|
{
|
||||||
|
if (strstr(update->integrate_style,"verlet"))
|
||||||
|
post_force(vflag);
|
||||||
|
else {
|
||||||
|
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
|
||||||
|
post_force_respa(vflag,nlevels_respa-1,0);
|
||||||
|
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void FixGLE::gle_integrate()
|
||||||
|
{
|
||||||
|
double **v = atom->v;
|
||||||
|
double *rmass = atom->rmass, smi, ismi;
|
||||||
|
int *type = atom->type;
|
||||||
|
int *mask = atom->mask;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||||
|
|
||||||
|
#ifdef GLE_DEBUG
|
||||||
|
printf("!MC! GLE THERMO STEP dt=%f\n",update->dt);
|
||||||
|
#endif
|
||||||
|
|
||||||
|
// loads momentum data (mass-scaled) into the temporary vectors for the propagation
|
||||||
|
int nk=0, ni=0; double deltae=0.0;
|
||||||
|
for (int i = 0; i < nlocal; i++) {
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
ni++;
|
||||||
|
if (rmass) {
|
||||||
|
smi = sqrt(rmass[i]);
|
||||||
|
} else {
|
||||||
|
smi=sqrt_m[type[i]];
|
||||||
|
}
|
||||||
|
for (int k = 0; k<3; k++) {
|
||||||
|
//also zeroes tmp1
|
||||||
|
gle_tmp1[nk]=0.0;
|
||||||
|
// first loads velocities and accumulates conserved quantity
|
||||||
|
gle_tmp2[nk]=v[i][k]*smi; deltae+= gle_tmp2[nk]*gle_tmp2[nk]; ++nk;
|
||||||
|
// then copies in the additional momenta
|
||||||
|
for (int j=0; j<ns; ++j)
|
||||||
|
gle_tmp2[nk++] = gle_s[i][k*ns+j];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// s(t+dt) = T.s ....
|
||||||
|
GLE::AkMult(ni*3,ns+1,ns+1,gle_tmp2,TT,gle_tmp1);
|
||||||
|
|
||||||
|
//fills up a vector of random numbers
|
||||||
|
for (int i = 0; i < 3*ni*(ns+1); i++) gle_tmp2[i] = random->gaussian();
|
||||||
|
|
||||||
|
// ... + S \xi
|
||||||
|
GLE::AkMult(ni*3,ns+1,ns+1,gle_tmp2,ST,gle_tmp1,1.0);
|
||||||
|
|
||||||
|
// unloads momentum data (mass-scaled) from the temporary vectors
|
||||||
|
nk=0;
|
||||||
|
for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) {
|
||||||
|
if (rmass) ismi = 1.0/sqrt(rmass[i]); else ismi=1.0/sqrt_m[type[i]];
|
||||||
|
for (int k = 0; k<3; k++)
|
||||||
|
{
|
||||||
|
// fetches new velocities and completes computation of the conserved quantity change
|
||||||
|
v[i][k]=gle_tmp1[nk]*ismi; deltae-= gle_tmp1[nk]*gle_tmp1[nk]; ++nk;
|
||||||
|
// stores the additional momenta in the gle_s buffer
|
||||||
|
for (int j=0; j<ns; ++j)
|
||||||
|
gle_s[i][k*ns+j]=gle_tmp1[nk++];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
energy += deltae*0.5*force->mvv2e;
|
||||||
|
}
|
||||||
|
|
||||||
|
void FixGLE::initial_integrate(int vflag)
|
||||||
|
{
|
||||||
|
double dtfm;
|
||||||
|
|
||||||
|
// update v and x of atoms in group
|
||||||
|
double **x = atom->x;
|
||||||
|
double **v = atom->v;
|
||||||
|
double **f = atom->f;
|
||||||
|
double *rmass = atom->rmass;
|
||||||
|
double *mass = atom->mass;
|
||||||
|
int *type = atom->type;
|
||||||
|
int *mask = atom->mask;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||||
|
|
||||||
|
gle_step--;
|
||||||
|
if (dogle && gle_step<1) gle_integrate();
|
||||||
|
|
||||||
|
#ifdef GLE_DEBUG
|
||||||
|
printf("!MC! GLE P1 STEP dt=%f\n",dtv);
|
||||||
|
printf("!MC! GLE Q STEP\n");
|
||||||
|
#endif
|
||||||
|
|
||||||
|
if (rmass) {
|
||||||
|
for (int i = 0; i < nlocal; i++)
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
dtfm = dtf / rmass[i];
|
||||||
|
v[i][0] += dtfm * f[i][0];
|
||||||
|
v[i][1] += dtfm * f[i][1];
|
||||||
|
v[i][2] += dtfm * f[i][2];
|
||||||
|
x[i][0] += dtv * v[i][0];
|
||||||
|
x[i][1] += dtv * v[i][1];
|
||||||
|
x[i][2] += dtv * v[i][2];
|
||||||
|
}
|
||||||
|
|
||||||
|
} else {
|
||||||
|
for (int i = 0; i < nlocal; i++)
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
dtfm = dtf / mass[type[i]];
|
||||||
|
v[i][0] += dtfm * f[i][0];
|
||||||
|
v[i][1] += dtfm * f[i][1];
|
||||||
|
v[i][2] += dtfm * f[i][2];
|
||||||
|
x[i][0] += dtv * v[i][0];
|
||||||
|
x[i][1] += dtv * v[i][1];
|
||||||
|
x[i][2] += dtv * v[i][2];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void FixGLE::final_integrate()
|
||||||
|
{
|
||||||
|
double dtfm;
|
||||||
|
|
||||||
|
// update v of atoms in group
|
||||||
|
|
||||||
|
double **v = atom->v;
|
||||||
|
double **f = atom->f;
|
||||||
|
double *rmass = atom->rmass;
|
||||||
|
double *mass = atom->mass;
|
||||||
|
int *type = atom->type;
|
||||||
|
int *mask = atom->mask;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||||
|
|
||||||
|
if (rmass) {
|
||||||
|
for (int i = 0; i < nlocal; i++)
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
dtfm = dtf / rmass[i];
|
||||||
|
v[i][0] += dtfm * f[i][0];
|
||||||
|
v[i][1] += dtfm * f[i][1];
|
||||||
|
v[i][2] += dtfm * f[i][2];
|
||||||
|
}
|
||||||
|
|
||||||
|
} else {
|
||||||
|
for (int i = 0; i < nlocal; i++)
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
dtfm = dtf / mass[type[i]];
|
||||||
|
v[i][0] += dtfm * f[i][0];
|
||||||
|
v[i][1] += dtfm * f[i][1];
|
||||||
|
v[i][2] += dtfm * f[i][2];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
#ifdef GLE_DEBUG
|
||||||
|
printf("!MC! GLE P2 STEP dt=%f\n",dtv);
|
||||||
|
#endif
|
||||||
|
if (dogle && gle_step<1) { gle_integrate(); gle_step=gle_every; }
|
||||||
|
|
||||||
|
// Change the temperature for the next step
|
||||||
|
double delta = update->ntimestep - update->beginstep;
|
||||||
|
delta /= update->endstep - update->beginstep;
|
||||||
|
t_target = t_start + delta * (t_stop - t_start);
|
||||||
|
if (t_stop != t_start && fnoneq == 0) {
|
||||||
|
// only updates if it is really necessary
|
||||||
|
const double kT = t_target * force->boltz / force->mvv2e;
|
||||||
|
memset(C,0,sizeof(double)*ns1sq);
|
||||||
|
for (int i=0; i<ns1sq; i+=(ns+2)) C[i]=kT;
|
||||||
|
init_gle();
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixGLE::initial_integrate_respa(int vflag, int ilevel, int iloop)
|
||||||
|
{
|
||||||
|
dtv = step_respa[ilevel];
|
||||||
|
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
|
||||||
|
|
||||||
|
// innermost level - NVE update of v and x
|
||||||
|
// all other levels - NVE update of v
|
||||||
|
|
||||||
|
if (ilevel==nlevels_respa-1) gle_integrate();
|
||||||
|
dogle=0;
|
||||||
|
if (ilevel == 0) initial_integrate(vflag);
|
||||||
|
else { final_integrate();}
|
||||||
|
}
|
||||||
|
|
||||||
|
void FixGLE::final_integrate_respa(int ilevel, int iloop)
|
||||||
|
{
|
||||||
|
dtv = step_respa[ilevel];
|
||||||
|
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
|
||||||
|
dogle=0;
|
||||||
|
final_integrate();
|
||||||
|
if (ilevel==nlevels_respa-1) gle_integrate();
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
double FixGLE::compute_scalar()
|
||||||
|
{
|
||||||
|
double energy_me = energy;
|
||||||
|
double energy_all;
|
||||||
|
MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||||
|
|
||||||
|
return energy_all;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
extract thermostat properties
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void *FixGLE::extract(const char *str, int &dim)
|
||||||
|
{
|
||||||
|
dim = 0;
|
||||||
|
if (strcmp(str,"t_target") == 0) {
|
||||||
|
return &t_target;
|
||||||
|
}
|
||||||
|
return NULL;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Called when a change to the target temperature is requested mid-run
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixGLE::reset_target(double t_new)
|
||||||
|
{
|
||||||
|
t_target = t_start = t_stop = t_new;
|
||||||
|
if (fnoneq == 0)
|
||||||
|
{
|
||||||
|
// only updates if it is really necessary
|
||||||
|
for (int i=0; i<ns1sq; ++i) C[i]=0.0;
|
||||||
|
for (int i=0; i<ns1sq; i+=(ns+2)) C[i]=t_target*force->boltz/force->mvv2e;
|
||||||
|
init_gle();
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
error->all(FLERR, "Cannot change temperature for a non-equilibrium GLE run");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Called when a change to the timestep is requested mid-run
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixGLE::reset_dt()
|
||||||
|
{
|
||||||
|
// set the time integration constants
|
||||||
|
dtv = update->dt;
|
||||||
|
dtf = 0.5 * update->dt * (force->ftm2v);
|
||||||
|
init_gle();
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
memory usage of local atom-based arrays
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
double FixGLE::memory_usage()
|
||||||
|
{
|
||||||
|
double bytes = atom->nmax*(3*ns+2*3*(ns+1))*sizeof(double);
|
||||||
|
return bytes;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
allocate local atom-based arrays
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixGLE::grow_arrays(int nmax)
|
||||||
|
{
|
||||||
|
memory->grow(gle_s, nmax, 3*ns,"gle:gle_s");
|
||||||
|
memory->grow(gle_tmp1, nmax*(ns+1)*3,"gle:tmp1");
|
||||||
|
memory->grow(gle_tmp2, nmax*(ns+1)*3,"gle:tmp2");
|
||||||
|
//zeroes out temporary buffers
|
||||||
|
for (int i=0; i< nmax*(ns+1)*3; ++i) gle_tmp1[i]=0.0;
|
||||||
|
for (int i=0; i< nmax*(ns+1)*3; ++i) gle_tmp2[i]=0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
copy values within local atom-based arrays
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixGLE::copy_arrays(int i, int j, int delflag)
|
||||||
|
{
|
||||||
|
for (int k = 0; k < 3*ns; k++) gle_s[j][k] = gle_s[i][k];
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Pack extended variables assoc. w/ atom i into buffer for exchange
|
||||||
|
with another processor
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixGLE::pack_exchange(int i, double *buf)
|
||||||
|
{
|
||||||
|
int m = 0;
|
||||||
|
for (int k = 0; k < 3*ns; k++) buf[m++] = gle_s[i][k];
|
||||||
|
return m;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Unpack extended variables from exchange with another processor
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixGLE::unpack_exchange(int nlocal, double *buf)
|
||||||
|
{
|
||||||
|
int m = 0;
|
||||||
|
for (int k = 0; k < 3*ns; k++) gle_s[nlocal][k] = buf[m++];
|
||||||
|
return m;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Pack extended variables assoc. w/ atom i into buffer for
|
||||||
|
writing to a restart file
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixGLE::pack_restart(int i, double *buf)
|
||||||
|
{
|
||||||
|
int m = 0;
|
||||||
|
buf[m++] = 3*ns + 1;
|
||||||
|
for (int k = 0; k < 3*ns; k=k+3)
|
||||||
|
{
|
||||||
|
buf[m++] = gle_s[i][k];
|
||||||
|
buf[m++] = gle_s[i][k+1];
|
||||||
|
buf[m++] = gle_s[i][k+2];
|
||||||
|
}
|
||||||
|
return m;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Unpack extended variables to restart the fix from a restart file
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixGLE::unpack_restart(int nlocal, int nth)
|
||||||
|
{
|
||||||
|
double **extra = atom->extra;
|
||||||
|
|
||||||
|
// skip to the nth set of extended variables
|
||||||
|
|
||||||
|
int m = 0;
|
||||||
|
for (int i = 0; i< nth; i++) m += static_cast<int> (extra[nlocal][m]);
|
||||||
|
m++;
|
||||||
|
|
||||||
|
for (int k = 0; k < 3*ns; k=k+3)
|
||||||
|
{
|
||||||
|
gle_s[nlocal][k] = extra[nlocal][m++];
|
||||||
|
gle_s[nlocal][k+1] = extra[nlocal][m++];
|
||||||
|
gle_s[nlocal][k+2] = extra[nlocal][m++];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Returns the number of items in atomic restart data associated with
|
||||||
|
local atom nlocal. Used in determining the total extra data stored by
|
||||||
|
fixes on a given processor.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixGLE::size_restart(int nlocal)
|
||||||
|
{
|
||||||
|
return 3*ns+1;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
Returns the maximum number of items in atomic restart data
|
||||||
|
Called in Modify::restart for peratom restart.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixGLE::maxsize_restart()
|
||||||
|
{
|
||||||
|
return 3*ns+1;
|
||||||
|
}
|
||||||
77
src/USER-MISC/fix_gle.h
Normal file
77
src/USER-MISC/fix_gle.h
Normal file
@ -0,0 +1,77 @@
|
|||||||
|
/* -*- c++ -*- ----------------------------------------------------------
|
||||||
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||||
|
http://lammps.sandia.gov, Sandia National Laboratories
|
||||||
|
Steve Plimpton, sjplimp@sandia.gov
|
||||||
|
|
||||||
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||||
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||||
|
certain rights in this software. This software is distributed under
|
||||||
|
the GNU General Public License.
|
||||||
|
|
||||||
|
See the README file in the top-level LAMMPS directory.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#ifdef FIX_CLASS
|
||||||
|
|
||||||
|
FixStyle(gle,FixGLE)
|
||||||
|
|
||||||
|
#else
|
||||||
|
|
||||||
|
#ifndef LMP_FIX_GLE_H
|
||||||
|
#define LMP_FIX_GLE_H
|
||||||
|
|
||||||
|
#include "fix.h"
|
||||||
|
|
||||||
|
namespace LAMMPS_NS {
|
||||||
|
|
||||||
|
class FixGLE : public Fix {
|
||||||
|
public:
|
||||||
|
FixGLE(class LAMMPS *, int, char **);
|
||||||
|
virtual ~FixGLE();
|
||||||
|
int setmask();
|
||||||
|
void init();
|
||||||
|
void setup(int);
|
||||||
|
void gle_integrate();
|
||||||
|
void initial_integrate_respa(int vflag, int ilevel, int iloop);
|
||||||
|
void final_integrate_respa(int ilevel, int iloop);
|
||||||
|
void initial_integrate(int vflag);
|
||||||
|
void final_integrate();
|
||||||
|
double compute_scalar();
|
||||||
|
void reset_target(double);
|
||||||
|
virtual void reset_dt();
|
||||||
|
|
||||||
|
double memory_usage();
|
||||||
|
void grow_arrays(int);
|
||||||
|
void copy_arrays(int, int, int);
|
||||||
|
int pack_exchange(int, double *);
|
||||||
|
int unpack_exchange(int, double *);
|
||||||
|
int pack_restart(int, double *);
|
||||||
|
void unpack_restart(int, int);
|
||||||
|
int size_restart(int);
|
||||||
|
int maxsize_restart();
|
||||||
|
|
||||||
|
virtual void *extract(const char *, int &);
|
||||||
|
|
||||||
|
void init_gle(); void init_gles();
|
||||||
|
protected:
|
||||||
|
int ns, ns1sq;
|
||||||
|
double *A, *C, *S, *T, *ST, *TT;
|
||||||
|
double *gle_tmp1, *gle_tmp2;
|
||||||
|
double t_start, t_stop, t_target;
|
||||||
|
double dtv, dtf;
|
||||||
|
|
||||||
|
int dogle, fnoneq, gle_every, gle_step;
|
||||||
|
class RanMars *random;
|
||||||
|
double *sqrt_m;
|
||||||
|
double *step_respa;
|
||||||
|
double energy;
|
||||||
|
int nlevels_respa;
|
||||||
|
|
||||||
|
double **gle_s;
|
||||||
|
double **vaux;
|
||||||
|
};
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#endif
|
||||||
437
src/USER-MISC/fix_ipi.cpp
Normal file
437
src/USER-MISC/fix_ipi.cpp
Normal file
@ -0,0 +1,437 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||||
|
http://lammps.sandia.gov, Sandia National Laboratories
|
||||||
|
Steve Plimpton, sjplimp@sandia.gov
|
||||||
|
|
||||||
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||||
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||||
|
certain rights in this software. This software is distributed under
|
||||||
|
the GNU General Public License.
|
||||||
|
|
||||||
|
See the README file in the top-level LAMMPS directory.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#include "mpi.h"
|
||||||
|
#include "stdio.h"
|
||||||
|
#include "string.h"
|
||||||
|
#include "fix_ipi.h"
|
||||||
|
#include "atom.h"
|
||||||
|
#include "force.h"
|
||||||
|
#include "update.h"
|
||||||
|
#include "respa.h"
|
||||||
|
#include "error.h"
|
||||||
|
#include "kspace.h"
|
||||||
|
#include "modify.h"
|
||||||
|
#include "compute.h"
|
||||||
|
#include "comm.h"
|
||||||
|
#include "neighbor.h"
|
||||||
|
#include "domain.h"
|
||||||
|
#include "compute_pressure.h"
|
||||||
|
#include "errno.h"
|
||||||
|
|
||||||
|
using namespace LAMMPS_NS;
|
||||||
|
using namespace FixConst;
|
||||||
|
|
||||||
|
/******************************************************************************************
|
||||||
|
* A fix to interface LAMMPS with i-PI - A Python interface for path integral molecular dynamics
|
||||||
|
* Michele Ceriotti, EPFL (2014)
|
||||||
|
* Please cite:
|
||||||
|
* Ceriotti, M., More, J., & Manolopoulos, D. E. (2014).
|
||||||
|
* i-PI: A Python interface for ab initio path integral molecular dynamics simulations.
|
||||||
|
* Computer Physics Communications, 185, 1019–1026. doi:10.1016/j.cpc.2013.10.027
|
||||||
|
* And see [http://github.com/i-pi/i-pi] to download a version of i-PI
|
||||||
|
******************************************************************************************/
|
||||||
|
|
||||||
|
// socket interface
|
||||||
|
#ifndef _WIN32
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <unistd.h>
|
||||||
|
#include <sys/types.h>
|
||||||
|
#include <sys/socket.h>
|
||||||
|
#include <netinet/in.h>
|
||||||
|
#include <sys/un.h>
|
||||||
|
#include <netdb.h>
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#define MSGLEN 12
|
||||||
|
|
||||||
|
/* Utility functions to simplify the interface with POSIX sockets */
|
||||||
|
|
||||||
|
static void open_socket(int &sockfd, int inet, int port, char* host,
|
||||||
|
Error *error)
|
||||||
|
/* Opens a socket.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
sockfd: The id of the socket that will be created.
|
||||||
|
inet: An integer that determines whether the socket will be an inet or unix
|
||||||
|
domain socket. Gives unix if 0, inet otherwise.
|
||||||
|
port: The port number for the socket to be created. Low numbers are often
|
||||||
|
reserved for important channels, so use of numbers of 4 or more digits is
|
||||||
|
recommended.
|
||||||
|
host: The name of the host server.
|
||||||
|
error: pointer to a LAMMPS Error object
|
||||||
|
*/
|
||||||
|
{
|
||||||
|
int ai_err;
|
||||||
|
|
||||||
|
#ifdef _WIN32
|
||||||
|
error->one(FLERR,"i-PI socket implementation requires UNIX environment");
|
||||||
|
#else
|
||||||
|
if (inet>0) { // creates an internet socket
|
||||||
|
|
||||||
|
// fetches information on the host
|
||||||
|
struct addrinfo hints, *res;
|
||||||
|
char service[256];
|
||||||
|
|
||||||
|
memset(&hints, 0, sizeof(hints));
|
||||||
|
hints.ai_socktype = SOCK_STREAM;
|
||||||
|
hints.ai_family = AF_UNSPEC;
|
||||||
|
hints.ai_flags = AI_PASSIVE;
|
||||||
|
|
||||||
|
sprintf(service,"%d",port); // convert the port number to a string
|
||||||
|
ai_err = getaddrinfo(host, service, &hints, &res);
|
||||||
|
if (ai_err!=0)
|
||||||
|
error->one(FLERR,"Error fetching host data. Wrong host name?");
|
||||||
|
|
||||||
|
// creates socket
|
||||||
|
sockfd = socket(res->ai_family, res->ai_socktype, res->ai_protocol);
|
||||||
|
if (sockfd < 0)
|
||||||
|
error->one(FLERR,"Error opening socket");
|
||||||
|
|
||||||
|
// makes connection
|
||||||
|
if (connect(sockfd, res->ai_addr, res->ai_addrlen) < 0)
|
||||||
|
error->one(FLERR,"Error opening INET socket: wrong port or server unreachable");
|
||||||
|
freeaddrinfo(res);
|
||||||
|
|
||||||
|
} else { // creates a unix socket
|
||||||
|
struct sockaddr_un serv_addr;
|
||||||
|
|
||||||
|
// fills up details of the socket addres
|
||||||
|
memset(&serv_addr, 0, sizeof(serv_addr));
|
||||||
|
serv_addr.sun_family = AF_UNIX;
|
||||||
|
strcpy(serv_addr.sun_path, "/tmp/ipi_");
|
||||||
|
strcpy(serv_addr.sun_path+9, host);
|
||||||
|
|
||||||
|
// creates the socket
|
||||||
|
sockfd = socket(AF_UNIX, SOCK_STREAM, 0);
|
||||||
|
|
||||||
|
// connects
|
||||||
|
if (connect(sockfd, (struct sockaddr *) &serv_addr, sizeof(serv_addr)) < 0)
|
||||||
|
error->one(FLERR,"Error opening UNIX socket: server may not be running "
|
||||||
|
"or the path to the socket unavailable");
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
}
|
||||||
|
|
||||||
|
static void writebuffer(int sockfd, const char *data, int len, Error* error)
|
||||||
|
/* Writes to a socket.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
sockfd: The id of the socket that will be written to.
|
||||||
|
data: The data to be written to the socket.
|
||||||
|
len: The length of the data in bytes.
|
||||||
|
*/
|
||||||
|
{
|
||||||
|
int n;
|
||||||
|
|
||||||
|
n = write(sockfd,data,len);
|
||||||
|
if (n < 0)
|
||||||
|
error->one(FLERR,"Error writing to socket: broken connection");
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
static void readbuffer(int sockfd, char *data, int len, Error* error)
|
||||||
|
/* Reads from a socket.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
sockfd: The id of the socket that will be read from.
|
||||||
|
data: The storage array for data read from the socket.
|
||||||
|
len: The length of the data in bytes.
|
||||||
|
*/
|
||||||
|
{
|
||||||
|
int n, nr;
|
||||||
|
|
||||||
|
n = nr = read(sockfd,data,len);
|
||||||
|
|
||||||
|
while (nr>0 && n<len ) {
|
||||||
|
nr=read(sockfd,&data[n],len-n);
|
||||||
|
n+=nr;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (n == 0)
|
||||||
|
error->one(FLERR,"Error reading from socket: broken connection");
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
FixIPI::FixIPI(LAMMPS *lmp, int narg, char **arg) :
|
||||||
|
Fix(lmp, narg, arg)
|
||||||
|
{
|
||||||
|
/* format for fix:
|
||||||
|
* fix num group_id ipi host port [unix]
|
||||||
|
*/
|
||||||
|
if (strcmp(style,"ipi") != 0 && narg < 5)
|
||||||
|
error->all(FLERR,"Illegal fix ipi command");
|
||||||
|
|
||||||
|
if (atom->tag_enable == 0)
|
||||||
|
error->all(FLERR,"Cannot use fix ipi without atom IDs");
|
||||||
|
|
||||||
|
if (atom->tag_consecutive() == 0)
|
||||||
|
error->all(FLERR,"Fix ipi requires consecutive atom IDs");
|
||||||
|
|
||||||
|
if (strcmp(arg[1],"all"))
|
||||||
|
error->warning(FLERR,"Fix ipi always uses group all");
|
||||||
|
|
||||||
|
strcpy(host, arg[3]);
|
||||||
|
port=force->inumeric(FLERR,arg[4]);
|
||||||
|
|
||||||
|
inet = ((narg > 5) && (strcmp(arg[5],"unix") ==0) ) ? 0 : 1;
|
||||||
|
master = (comm->me==0) ? 1 : 0;
|
||||||
|
hasdata = bsize = 0;
|
||||||
|
|
||||||
|
// creates a temperature compute for all atoms
|
||||||
|
char** newarg = new char*[3];
|
||||||
|
newarg[0] = (char *) "IPI_TEMP";
|
||||||
|
newarg[1] = (char *) "all";
|
||||||
|
newarg[2] = (char *) "temp";
|
||||||
|
modify->add_compute(3,newarg);
|
||||||
|
delete [] newarg;
|
||||||
|
|
||||||
|
// creates a pressure compute to extract the virial
|
||||||
|
newarg = new char*[5];
|
||||||
|
newarg[0] = (char *) "IPI_PRESS";
|
||||||
|
newarg[1] = (char *) "all";
|
||||||
|
newarg[2] = (char *) "pressure";
|
||||||
|
newarg[3] = (char *) "IPI_TEMP";
|
||||||
|
newarg[4] = (char *) "virial";
|
||||||
|
modify->add_compute(5,newarg);
|
||||||
|
delete [] newarg;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
FixIPI::~FixIPI()
|
||||||
|
{
|
||||||
|
if (bsize) delete[] buffer;
|
||||||
|
modify->delete_compute("IPI_TEMP");
|
||||||
|
modify->delete_compute("IPI_PRESS");
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int FixIPI::setmask()
|
||||||
|
{
|
||||||
|
int mask = 0;
|
||||||
|
mask |= INITIAL_INTEGRATE;
|
||||||
|
mask |= FINAL_INTEGRATE;
|
||||||
|
return mask;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void FixIPI::init()
|
||||||
|
{
|
||||||
|
//only opens socket on master process
|
||||||
|
if (master) open_socket(ipisock, inet, port, host, error);
|
||||||
|
else ipisock=0;
|
||||||
|
//! should check for success in socket opening -- but the current open_socket routine dies brutally if unsuccessful
|
||||||
|
|
||||||
|
// asks for evaluation of PE at first step
|
||||||
|
modify->compute[modify->find_compute("thermo_pe")]->invoked_scalar = -1;
|
||||||
|
modify->addstep_compute_all(update->ntimestep + 1);
|
||||||
|
|
||||||
|
kspace_flag = (force->kspace) ? 1 : 0;
|
||||||
|
|
||||||
|
// makes sure that neighbor lists are re-built at each step (cannot make assumptions when cycling over beads!)
|
||||||
|
neighbor->delay = 0;
|
||||||
|
neighbor->every = 1;
|
||||||
|
}
|
||||||
|
|
||||||
|
void FixIPI::initial_integrate(int vflag)
|
||||||
|
{
|
||||||
|
/* This is called at the beginning of the integration loop,
|
||||||
|
* and will be used to read positions from the socket. Then,
|
||||||
|
* everything should be updated, since there is no guarantee
|
||||||
|
* that successive snapshots will be close together (think
|
||||||
|
* of parallel tempering for instance) */
|
||||||
|
|
||||||
|
char header[MSGLEN+1];
|
||||||
|
|
||||||
|
if (hasdata)
|
||||||
|
error->all(FLERR, "i-PI got out of sync in initial_integrate and will die!");
|
||||||
|
|
||||||
|
double cellh[9], cellih[9];
|
||||||
|
int nat;
|
||||||
|
if (master) { // only read positions on master
|
||||||
|
|
||||||
|
// wait until something happens
|
||||||
|
while (true) {
|
||||||
|
// while i-PI just asks for status, signal we are ready and wait
|
||||||
|
readbuffer(ipisock, header, MSGLEN, error); header[MSGLEN]=0;
|
||||||
|
|
||||||
|
if (strcmp(header,"STATUS ") == 0 )
|
||||||
|
writebuffer(ipisock,"READY ",MSGLEN, error);
|
||||||
|
else break;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (strcmp(header,"EXIT ") == 0 )
|
||||||
|
error->one(FLERR, "Got EXIT message from i-PI. Now leaving!");
|
||||||
|
|
||||||
|
// when i-PI signals it has positions to evaluate new forces,
|
||||||
|
// read positions and cell data
|
||||||
|
if (strcmp(header,"POSDATA ") == 0 ) {
|
||||||
|
readbuffer(ipisock, (char*) cellh, 9*8, error);
|
||||||
|
readbuffer(ipisock, (char*) cellih, 9*8, error);
|
||||||
|
readbuffer(ipisock, (char*) &nat, 4, error);
|
||||||
|
|
||||||
|
// allocate buffer, but only do this once.
|
||||||
|
if (bsize==0) {
|
||||||
|
bsize=3*nat;
|
||||||
|
buffer = new double[bsize];
|
||||||
|
} else if (bsize != 3*nat)
|
||||||
|
error->one(FLERR, "Number of atoms changed along the way.");
|
||||||
|
|
||||||
|
// finally read position data into buffer
|
||||||
|
readbuffer(ipisock, (char*) buffer, 8*bsize, error);
|
||||||
|
} else
|
||||||
|
error->one(FLERR, "Wrapper did not send positions, I will now die!");
|
||||||
|
}
|
||||||
|
|
||||||
|
// shares the atomic coordinates with everyone
|
||||||
|
MPI_Bcast(&nat,1,MPI_INT,0,world);
|
||||||
|
// must also allocate the buffer on the non-head nodes
|
||||||
|
if (bsize==0) {
|
||||||
|
bsize=3*nat;
|
||||||
|
buffer = new double[bsize];
|
||||||
|
}
|
||||||
|
MPI_Bcast(cellh,9,MPI_DOUBLE,0,world);
|
||||||
|
MPI_Bcast(cellih,9,MPI_DOUBLE,0,world);
|
||||||
|
MPI_Bcast(buffer,bsize,MPI_DOUBLE,0,world);
|
||||||
|
|
||||||
|
//updates atomic coordinates and cell based on the data received
|
||||||
|
double *boxhi = domain->boxhi;
|
||||||
|
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;
|
||||||
|
domain->xy = cellh[1]*posconv;
|
||||||
|
domain->xz = cellh[2]*posconv;
|
||||||
|
domain->yz = cellh[5]*posconv;
|
||||||
|
|
||||||
|
// signal that the box has (or may have) changed
|
||||||
|
domain->reset_box();
|
||||||
|
domain->box_change = 1;
|
||||||
|
if (kspace_flag) force->kspace->setup();
|
||||||
|
|
||||||
|
// picks local atoms from the buffer
|
||||||
|
double **x = atom->x;
|
||||||
|
int *mask = atom->mask;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||||
|
for (int i = 0; i < nlocal; i++) {
|
||||||
|
if (mask[i] & groupbit) {
|
||||||
|
x[i][0]=buffer[3*(atom->tag[i]-1)+0]*posconv;
|
||||||
|
x[i][1]=buffer[3*(atom->tag[i]-1)+1]*posconv;
|
||||||
|
x[i][2]=buffer[3*(atom->tag[i]-1)+2]*posconv;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// compute PE. makes sure that it will be evaluated at next step
|
||||||
|
modify->compute[modify->find_compute("thermo_pe")]->invoked_scalar = -1;
|
||||||
|
modify->addstep_compute_all(update->ntimestep+1);
|
||||||
|
|
||||||
|
hasdata=1;
|
||||||
|
}
|
||||||
|
|
||||||
|
void FixIPI::final_integrate()
|
||||||
|
{
|
||||||
|
/* This is called after forces and energy have been computed. Now we only need to
|
||||||
|
* communicate them back to i-PI so that the integration can continue. */
|
||||||
|
char header[MSGLEN+1];
|
||||||
|
double vir[9], pot=0.0;
|
||||||
|
double forceconv, potconv, posconv, pressconv, posconv3;
|
||||||
|
char retstr[1024];
|
||||||
|
|
||||||
|
// conversions from LAMMPS units to atomic units, which are used by i-PI
|
||||||
|
potconv=3.1668152e-06/force->boltz;
|
||||||
|
posconv=0.52917721*force->angstrom;
|
||||||
|
posconv3=posconv*posconv*posconv;
|
||||||
|
forceconv=potconv*posconv;
|
||||||
|
pressconv=1/force->nktv2p*potconv*posconv3;
|
||||||
|
|
||||||
|
// compute for potential energy
|
||||||
|
pot=modify->compute[modify->find_compute("thermo_pe")]->compute_scalar();
|
||||||
|
pot*=potconv;
|
||||||
|
|
||||||
|
// probably useless check
|
||||||
|
if (!hasdata)
|
||||||
|
error->all(FLERR, "i-PI got out of sync in final_integrate and will die!");
|
||||||
|
|
||||||
|
int nat=bsize/3;
|
||||||
|
double **f= atom->f;
|
||||||
|
double lbuf[bsize];
|
||||||
|
|
||||||
|
// reassembles the force vector from the local arrays
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||||
|
for (int i = 0; i < bsize; ++i) lbuf[i]=0.0;
|
||||||
|
for (int i = 0; i < nlocal; i++) {
|
||||||
|
lbuf[3*(atom->tag[i]-1)+0]=f[i][0]*forceconv;
|
||||||
|
lbuf[3*(atom->tag[i]-1)+1]=f[i][1]*forceconv;
|
||||||
|
lbuf[3*(atom->tag[i]-1)+2]=f[i][2]*forceconv;
|
||||||
|
}
|
||||||
|
MPI_Allreduce(lbuf,buffer,bsize,MPI_DOUBLE,MPI_SUM,world);
|
||||||
|
|
||||||
|
for (int i = 0; i < 9; ++i) vir[i]=0.0;
|
||||||
|
|
||||||
|
int press_id = modify->find_compute("IPI_PRESS");
|
||||||
|
Compute* comp_p = modify->compute[press_id];
|
||||||
|
comp_p->compute_vector();
|
||||||
|
double myvol = domain->xprd*domain->yprd*domain->zprd/posconv3;
|
||||||
|
|
||||||
|
vir[0] = comp_p->vector[0]*pressconv*myvol;
|
||||||
|
vir[4] = comp_p->vector[1]*pressconv*myvol;
|
||||||
|
vir[8] = comp_p->vector[2]*pressconv*myvol;
|
||||||
|
vir[1] = comp_p->vector[3]*pressconv*myvol;
|
||||||
|
vir[2] = comp_p->vector[4]*pressconv*myvol;
|
||||||
|
vir[5] = comp_p->vector[5]*pressconv*myvol;
|
||||||
|
retstr[0]=0;
|
||||||
|
|
||||||
|
if (master) {
|
||||||
|
while (true) {
|
||||||
|
readbuffer(ipisock, header, MSGLEN, error); header[MSGLEN]=0;
|
||||||
|
|
||||||
|
if (strcmp(header,"STATUS ") == 0 )
|
||||||
|
writebuffer(ipisock,"HAVEDATA ",MSGLEN, error);
|
||||||
|
else break;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (strcmp(header,"EXIT ") == 0 )
|
||||||
|
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);
|
||||||
|
writebuffer(ipisock,(char*) buffer, bsize*8, error);
|
||||||
|
writebuffer(ipisock,(char*) vir,9*8, error);
|
||||||
|
nat=strlen(retstr); writebuffer(ipisock,(char*) &nat,4, error);
|
||||||
|
writebuffer(ipisock,(char*) retstr, nat, error);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
error->one(FLERR, "Wrapper did not ask for forces, I will now die!");
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
hasdata=0;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
45
src/USER-MISC/fix_ipi.h
Normal file
45
src/USER-MISC/fix_ipi.h
Normal file
@ -0,0 +1,45 @@
|
|||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||||
|
http://lammps.sandia.gov, Sandia National Laboratories
|
||||||
|
Steve Plimpton, sjplimp@sandia.gov
|
||||||
|
|
||||||
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||||
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||||
|
certain rights in this software. This software is distributed under
|
||||||
|
the GNU General Public License.
|
||||||
|
|
||||||
|
See the README file in the top-level LAMMPS directory.
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#ifdef FIX_CLASS
|
||||||
|
|
||||||
|
FixStyle(ipi,FixIPI)
|
||||||
|
|
||||||
|
#else
|
||||||
|
|
||||||
|
#ifndef LMP_FIX_IPI_H
|
||||||
|
#define LMP_FIX_IPI_H
|
||||||
|
|
||||||
|
#include "fix.h"
|
||||||
|
|
||||||
|
namespace LAMMPS_NS {
|
||||||
|
|
||||||
|
class FixIPI : public Fix {
|
||||||
|
public:
|
||||||
|
FixIPI(class LAMMPS *, int, char **);
|
||||||
|
virtual ~FixIPI();
|
||||||
|
int setmask();
|
||||||
|
virtual void init();
|
||||||
|
virtual void initial_integrate(int);
|
||||||
|
virtual void final_integrate();
|
||||||
|
|
||||||
|
protected:
|
||||||
|
char host[1024]; int port; int inet, master, hasdata;
|
||||||
|
int ipisock, me; double *buffer; long bsize;
|
||||||
|
int kspace_flag;
|
||||||
|
};
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#endif
|
||||||
Reference in New Issue
Block a user