multi-integrator version of adaptglok: semi-implicit Euler, Verlet and Leap Frog (89c9c32c) · Commits · 郑智淋 / lammps

src/min.cpp

+8 −0

Changes for src/min.cpp: 8 added lines, 0 removed lines.
Original line number	Diff line number	Diff line
		@@ -62,6 +62,7 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
		alpha_shrink = 0.99;
		tmax = 2.0;
		tmin = 0.02;
		integrator = 0;

		elist_global = elist_atom = NULL;
		vlist_global = vlist_atom = NULL;
		@@ -679,6 +680,13 @@ void Min::modify_params(int narg, char **arg)
		if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
		tmin = force->numeric(FLERR,arg[iarg+1]);
		iarg += 2;
		} else if (strcmp(arg[iarg],"integrator") == 0) {
		if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
		if (strcmp(arg[iarg+1],"euler") == 0) integrator = 0;
		else if (strcmp(arg[iarg+1],"verlet") == 0) integrator = 1;
		else if (strcmp(arg[iarg+1],"leapfrog") == 0) integrator = 2;
		else error->all(FLERR,"Illegal min_modify command");
		iarg += 2;
		} else if (strcmp(arg[iarg],"line") == 0) {
		if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
		if (strcmp(arg[iarg+1],"backtrack") == 0) linestyle = 0;

src/min.h

+1 −0

Changes for src/min.h: 1 added line, 0 removed lines.
Original line number	Diff line number	Diff line
		@@ -62,6 +62,7 @@ class Min : protected Pointers {
		double dt_grow,dt_shrink; // timestep increase, decrease (adaptglok)
		double alpha0,alpha_shrink; // mixing velocities+forces coefficient (adaptglok)
		double tmax,tmin; // timestep max, min (adaptglok)
		int integrator; // choose the style of time integrator (adaptglok)

		int nelist_global,nelist_atom; // # of PE,virial computes to check
		int nvlist_global,nvlist_atom;

src/min_adaptglok.cpp

+120 −5

Changes for src/min_adaptglok.cpp: 120 added lines, 5 removed lines.
Original line number	Diff line number	Diff line
		@@ -79,14 +79,37 @@ int MinAdaptGlok::iterate(int maxiter)
		double dtvone,dtv,dtf,dtfm;
		int flag,flagall;

		double *rmass = atom->rmass;
		double *mass = atom->mass;
		int *type = atom->type;
		alpha_final = 0.0;

		double **f = atom->f;
		double **v = atom->v;

		alpha_final = 0.0;
		// Leap Frog integration initialization

		if (integrator == 2) {

		double *rmass = atom->rmass;
		double *mass = atom->mass;
		int *type = atom->type;
		int nlocal = atom->nlocal;
		dtf = 0.5 * dt * force->ftm2v;
		if (rmass) {
		for (int i = 0; i < nlocal; i++) {
		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++) {
		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];
		}
		}

		}

		for (int iter = 0; iter < maxiter; iter++) {

		@@ -98,6 +121,12 @@ int MinAdaptGlok::iterate(int maxiter)

		// vdotfall = v dot f

		// Euler \|\| Leap Frog integration

		if (integrator == 0 \|\| integrator == 2) {
		double **v = atom->v;
		double **f = atom->f;
		}
		int nlocal = atom->nlocal;

		vdotf = 0.0;
		@@ -187,6 +216,10 @@ int MinAdaptGlok::iterate(int maxiter)

		// limit timestep so no particle moves further than dmax

		double *rmass = atom->rmass;
		double *mass = atom->mass;
		int *type = atom->type;

		dtvone = dt;

		for (int i = 0; i < nlocal; i++) {
		@@ -204,9 +237,51 @@ int MinAdaptGlok::iterate(int maxiter)
		MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld);
		}

		double **x = atom->x;

		// Semi-implicit Euler integration

		if (integrator == 0) {

		dtf = dtv * force->ftm2v;

		if (rmass) {
		for (int i = 0; i < nlocal; i++) {
		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];
		if (vdotfall > 0.0) {
		v[i][0] = scale1v[i][0] + scale2f[i][0];
		v[i][1] = scale1v[i][1] + scale2f[i][1];
		v[i][2] = scale1v[i][2] + scale2f[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++) {
		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];
		if (vdotfall > 0.0) {
		v[i][0] = scale1v[i][0] + scale2f[i][0];
		v[i][1] = scale1v[i][1] + scale2f[i][1];
		v[i][2] = scale1v[i][2] + scale2f[i][2];
		}
		x[i][0] += dtv * v[i][0];
		x[i][1] += dtv * v[i][1];
		x[i][2] += dtv * v[i][2];
		}
		}


		// Velocity Verlet integration

		double **x = atom->x;
		}else if (integrator == 1) {

		dtf = 0.5 * dtv * force->ftm2v;

		if (rmass) {
		@@ -260,6 +335,46 @@ int MinAdaptGlok::iterate(int maxiter)
		}
		}

		// Leap Frog integration

		}else if (integrator == 2) {

		dtf = dtv * force->ftm2v;

		if (rmass) {
		for (int i = 0; i < nlocal; i++) {
		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];
		if (vdotfall > 0.0) {
		v[i][0] = scale1v[i][0] + scale2f[i][0];
		v[i][1] = scale1v[i][1] + scale2f[i][1];
		v[i][2] = scale1v[i][2] + scale2f[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++) {
		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];
		if (vdotfall > 0.0) {
		v[i][0] = scale1v[i][0] + scale2f[i][0];
		v[i][1] = scale1v[i][1] + scale2f[i][1];
		v[i][2] = scale1v[i][2] + scale2f[i][2];
		}
		x[i][0] += dtv * v[i][0];
		x[i][1] += dtv * v[i][1];
		x[i][2] += dtv * v[i][2];
		}
		}

		}

		eprevious = ecurrent;
		ecurrent = energy_force(0);
		neval++;

Admin message