jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added potential gradient columns 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5117 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps 2010-10-23 00:05:44 +00:00
parent 3d4e086b2d
commit 86422df675
3 changed files with 53 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/REPLICA/fix_neb.cpp
View File

@ -168,6 +168,17 @@ void FixNEB::min_post_force(int vflag)
pe->addstep(update->ntimestep+1);
// Compute norm of GradV for log output
double **f = atom->f;
double fsq = 0.0;
for (int i = 0; i < nlocal; i++) {
fsq += f[i][0]*f[i][0]+f[i][1]*f[i][1]+f[i][2]*f[i][2];
}
MPI_Allreduce(&fsq,&gradvnorm,1,MPI_DOUBLE,MPI_MAX,world);
gradvnorm = sqrt(gradvnorm);
// if this is first or last replica, no change to forces, just return
if (ireplica == 0 || ireplica == nreplica-1) {
@ -274,13 +285,12 @@ void FixNEB::min_post_force(int vflag)
// see Henkelman & Jonsson 2000 paper, eqs 3,4,12
// see Henkelman & Jonsson 2000a paper, eq 5
double **f = atom->f;
double dot = 0.0;
for (int i = 0; i < nlocal; i++)
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] +
f[i][2]*tangent[i][2];
}
double prefactor;
if (ireplica == rclimber) prefactor = -2.0*dot;

1
src/REPLICA/fix_neb.h
View File

@ -28,6 +28,7 @@ class FixNEB : public Fix {
public:
double veng,plen,nlen;
int rclimber;
double gradvnorm;
FixNEB(class LAMMPS *, int, char **);
~FixNEB();

40
src/REPLICA/neb.cpp
View File

@ -138,9 +138,11 @@ void NEB::command(int narg, char **arg)
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
}
print_status();
@ -203,9 +205,11 @@ void NEB::command(int narg, char **arg)
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
}
print_status();
@ -381,11 +385,43 @@ void NEB::print_status()
double endpt = rdist[nreplica-1] = rdist[nreplica-2] + all[nreplica-2][2];
for (int i = 1; i < nreplica; i++)
rdist[i] /= endpt;
// look up GradV for the initial, final, and climbing replicas
// These are identical to fnorm2, but better to be safe we
// take them straight from fix_neb
double gradvnorm0, gradvnorm1, gradvnormc;
int irep;
irep = 0;
if (me_universe == irep) gradvnorm0 = fneb->gradvnorm;
MPI_Bcast(&gradvnorm0,1,MPI_DOUBLE,irep,uworld);
irep = nreplica-1;
if (me_universe == irep) gradvnorm1 = fneb->gradvnorm;
MPI_Bcast(&gradvnorm1,1,MPI_DOUBLE,irep,uworld);
irep = fneb->rclimber;
if (irep > -1) {
if (me_universe == irep) gradvnormc = fneb->gradvnorm;
MPI_Bcast(&gradvnormc,1,MPI_DOUBLE,irep,uworld);
} else {
double vmax = all[0][0];
int top = 0;
for (int m = 1; m < nreplica; m++)
if (vmax < all[m][0]) {
vmax = all[m][0];
top = m;
}
irep = top;
if (me_universe == irep) gradvnormc = fneb->gradvnorm;
MPI_Bcast(&gradvnormc,1,MPI_DOUBLE,irep,uworld);
}
if (me_universe == 0) {
if (universe->uscreen) {
fprintf(universe->uscreen,"%d %g %g ",update->ntimestep,
fmaxreplica,fmaxatom);
fprintf(universe->uscreen,"%g %g %g ",
gradvnorm0,gradvnorm1,gradvnormc);
for (int i = 0; i < nreplica; i++)
fprintf(universe->uscreen,"%g %g ",rdist[i],all[i][0]);
fprintf(universe->uscreen,"\n");
@ -393,6 +429,8 @@ void NEB::print_status()
if (universe->ulogfile) {
fprintf(universe->ulogfile,"%d %g %g ",update->ntimestep,
fmaxreplica,fmaxatom);
fprintf(universe->ulogfile,"%g %g %g ",
gradvnorm0,gradvnorm1,gradvnormc);
for (int i = 0; i < nreplica; i++)
fprintf(universe->ulogfile,"%g %g ",rdist[i],all[i][0]);
fprintf(universe->ulogfile,"\n");