ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Commit JT 031319

Browse Source 
- improved gneb functions
- correct name in fix_neb (Weinan)
This commit is contained in:
julient31
2019-03-13 17:02:40 -06:00
parent f6fb8b220d
commit 8c50c3d7c8
5 changed files with 82 additions and 420 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
doc/src/fix_neb.txt
View File

@ -97,7 +97,7 @@ Note that in this case the specified {Kspring} is in force/distance
units.
With a value of {ideal}, the spring force is computed as suggested in
"(WeinenE)"_#WeinenE :
"(WeinanE)"_#WeinanE :
Fnudge_parallel = -{Kspring} * (RD-RDideal) / (2 * meanDist) :pre
@ -224,8 +224,8 @@ specified (no inter-replica force on the end replicas).
[(Henkelman2)] Henkelman, Uberuaga, Jonsson, J Chem Phys, 113,
9901-9904 (2000).
:link(WeinenE)
[(WeinenE)] E, Ren, Vanden-Eijnden, Phys Rev B, 66, 052301 (2002).
:link(WeinanE)
[(WeinanE)] E, Ren, Vanden-Eijnden, Phys Rev B, 66, 052301 (2002).
:link(Jonsson)
[(Jonsson)] Jonsson, Mills and Jacobsen, in Classical and Quantum

1
doc/src/min_modify.txt
View File

@ -71,6 +71,7 @@ could move in the gradient direction to reduce forces further.
Keywords {alpha_damp} and {discret_factor} only make sense when
a {spinmin} minimization style is declared.
Default values are alpha_damp = 1.0 and discret_factor = 10.0.
[Restrictions:] none

72
src/REPLICA/fix_neb_spin.cpp
View File

@ -69,65 +69,25 @@ FixNEB_spin::FixNEB_spin(LAMMPS *lmp, int narg, char **arg) :
kspringPerp = 0.0;
kspringIni = 1.0;
kspringFinal = 1.0;
// only regular neb for now
SpinLattice = false;
SpinLattice = false; // no spin-lattice neb for now
// no available fix neb/spin options for now
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"lattice") == 0)
error->all(FLERR,"Illegal fix neb command");
}
/*
if (strcmp(arg[iarg],"parallel") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix neb command");
if (strcmp(arg[iarg+1],"ideal") == 0) {
NEBLongRange = true;
StandardNEB = false;
} else if (strcmp(arg[iarg+1],"neigh") == 0) {
NEBLongRange = false;
StandardNEB = true;
} else error->all(FLERR,"Illegal fix neb command");
error->all(FLERR,"Illegal fix neb command");
iarg += 2;
} else if (strcmp(arg[iarg],"perp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix neb command");
PerpSpring = true;
kspringPerp = force->numeric(FLERR,arg[iarg+1]);
if (kspringPerp == 0.0) PerpSpring = false;
if (kspringPerp < 0.0) error->all(FLERR,"Illegal fix neb command");
error->all(FLERR,"Illegal fix neb command");
iarg += 2;
} else if (strcmp (arg[iarg],"end") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix neb command");
if (strcmp(arg[iarg+1],"first") == 0) {
FreeEndIni = true;
kspringIni = force->numeric(FLERR,arg[iarg+2]);
} else if (strcmp(arg[iarg+1],"last") == 0) {
FreeEndFinal = true;
FinalAndInterWithRespToEIni = false;
FreeEndFinalWithRespToEIni = false;
kspringFinal = force->numeric(FLERR,arg[iarg+2]);
} else if (strcmp(arg[iarg+1],"last/efirst") == 0) {
FreeEndFinal = false;
FinalAndInterWithRespToEIni = false;
FreeEndFinalWithRespToEIni = true;
kspringFinal = force->numeric(FLERR,arg[iarg+2]);
} else if (strcmp(arg[iarg+1],"last/efirst/middle") == 0) {
FreeEndFinal = false;
FinalAndInterWithRespToEIni = true;
FreeEndFinalWithRespToEIni = true;
kspringFinal = force->numeric(FLERR,arg[iarg+2]);
} else error->all(FLERR,"Illegal fix neb command");
iarg += 3;
} else if (strcmp (arg[iarg],"lattice") == 0) {
iarg += 2;
} else error->all(FLERR,"Illegal fix neb command");
}
*/
// nreplica = number of partitions
// ireplica = which world I am in universe
// nprocs_universe = # of procs in all replicase
@ -172,10 +132,6 @@ FixNEB_spin::FixNEB_spin(LAMMPS *lmp, int narg, char **arg) :
modify->add_compute(3,newarg);
delete [] newarg;
// might need a test
// => check if pe does not compute mech potentials
// initialize local storage
maxlocal = -1;
@ -374,14 +330,12 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
pe->addstep(update->ntimestep+1);
double **x = atom->x;
// spin quantities
double **sp = atom->sp;
int *mask = atom->mask;
double dot = 0.0;
double prefactor = 0.0;
double **f = atom->f;
// spin quantities
//double **f = atom->f;
double **fm = atom->fm;
int nlocal = atom->nlocal;
@ -397,9 +351,10 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
// computation of the tangent vector
// final replica
if (ireplica == nreplica-1) {
// final replica
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@ -470,8 +425,10 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
//}
}
// initial replica
} else if (ireplica == 0) {
// initial replica
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
@ -537,7 +494,6 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
//}
}
// in-between replica
} else {
// not the first or last replica
@ -770,7 +726,7 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
//prefactor = -2.0*dot;
} else {
if (NEBLongRange) {
error->all(FLERR,"NEB_spin climber option not yet active");
error->all(FLERR,"Long Range NEB_spin climber option not yet active");
//prefactor = -dot - kspring*(lenuntilIm-idealPos)/(2*meanDist);
} else if (StandardNEB) {
prefactor = -dot + kspring*(nlen-plen);

417
src/REPLICA/neb_spin.cpp
View File

@ -72,7 +72,8 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
double delspx,delspy,delspz;
etol = etol_in;
ftol = ftol_in;
//ftol = ftol_in;
ttol = ftol_in;
n1steps = n1steps_in;
n2steps = n2steps_in;
nevery = nevery_in;
@ -85,12 +86,10 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
uworld = universe->uworld;
MPI_Comm_rank(world,&me);
// generate linear interpolate replica
// fraction to interpolate intermediate replica
double fraction = ireplica/(nreplica-1.0);
double **x = atom->x;
// spin quantitites
double **sp = atom->sp;
int nlocal = atom->nlocal;
@ -111,19 +110,11 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
sp[i][1] = spfinal[1];
sp[i][2] = spfinal[2];
//delx = buf_final[ii] - buf_init[ii];
//dely = buf_final[ii+1] - buf_init[ii+1];
//delz = buf_final[ii+2] - buf_init[ii+2];
// adjust distance if pbc
// not implemented yet
//domain->minimum_image(delx,dely,delz);
// need to define a procedure for circular initialization
//x[i][0] = buf_init[ii] + fraction*delx;
//x[i][1] = buf_init[ii+1] + fraction*dely;
//x[i][2] = buf_init[ii+2] + fraction*delz;
// need to define a better procedure for circular initialization
ii += 3;
}
}
@ -143,22 +134,14 @@ NEB_spin::~NEB_spin()
void NEB_spin::command(int narg, char **arg)
{
//printf("test 1 \n");
// test 1
double **sp1 = atom->sp;
//printf("test 1 atom: i=%d,%g,%g,%g \n",1,sp1[1][0],sp1[1][1],sp1[1][2]);
//error->all(FLERR,"end neb_spin test");
if (domain->box_exist == 0)
error->all(FLERR,"NEB_spin command before simulation box is defined");
if (narg < 6) error->universe_all(FLERR,"Illegal NEB_spin command");
etol = force->numeric(FLERR,arg[0]);
ftol = force->numeric(FLERR,arg[1]);
//ftol = force->numeric(FLERR,arg[1]);
ttol = force->numeric(FLERR,arg[1]);
n1steps = force->inumeric(FLERR,arg[2]);
n2steps = force->inumeric(FLERR,arg[3]);
nevery = force->inumeric(FLERR,arg[4]);
@ -166,7 +149,8 @@ void NEB_spin::command(int narg, char **arg)
// error checks
if (etol < 0.0) error->all(FLERR,"Illegal NEB_spin command");
if (ftol < 0.0) error->all(FLERR,"Illegal NEB_spin command");
//if (ftol < 0.0) error->all(FLERR,"Illegal NEB_spin command");
if (ttol < 0.0) error->all(FLERR,"Illegal NEB_spin command");
if (nevery <= 0) error->universe_all(FLERR,"Illegal NEB_spin command");
if (n1steps % nevery || n2steps % nevery)
error->universe_all(FLERR,"Illegal NEB_spin command");
@ -209,14 +193,6 @@ void NEB_spin::command(int narg, char **arg)
verbose=false;
if (strcmp(arg[narg-1],"verbose") == 0) verbose=true;
// test 1
double **sp = atom->sp;
//printf("test 2 atom: i=%d,%g,%g,%g \n",1,sp[1][0],sp[1][1],sp[1][2]);
//error->all(FLERR,"end neb_spin test");
// run the NEB_spin calculation
run();
}
@ -250,22 +226,24 @@ void NEB_spin::run()
update->whichflag = 2;
update->etol = etol;
update->ftol = ftol;
//update->ftol = ftol;
update->ftol = ttol; // update->ftol is a torque tolerance
update->multireplica = 1;
lmp->init();
// put flag to check gilbert damping procedure is set
// check if correct minimizer is setup
//if (update->minimize->searchflag)
// error->all(FLERR,"NEB_spin requires damped dynamics minimizer");
if (update->minimize->searchflag)
error->all(FLERR,"NEB_spin requires damped dynamics minimizer");
if (strcmp(update->minimize_style,"spinmin") != 0)
error->all(FLERR,"NEB_spin requires spinmin minimizer");
// setup regular NEB_spin minimization
FILE *uscreen = universe->uscreen;
FILE *ulogfile = universe->ulogfile;
//printf("test before run 1 \n");
if (me_universe == 0 && uscreen)
fprintf(uscreen,"Setting up regular NEB_spin ...\n");
@ -276,22 +254,23 @@ void NEB_spin::run()
if (update->laststep < 0)
error->all(FLERR,"Too many timesteps for NEB_spin");
//printf("test before run 2 \n");
update->minimize->setup();
//printf("test before run 3 \n");
if (me_universe == 0) {
if (uscreen) {
if (verbose) {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
fprintf(uscreen,"Step MaxReplicaTorque MaxAtomTorque "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
"... ReplicaTorqueN MaxAtomTorqueN\n");
//fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
// "GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
// "RDN PEN pathangle1 angletangrad1 anglegrad1 gradV1 "
// "ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
// "... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
fprintf(uscreen,"Step MaxReplicaTorque MaxAtomTorque "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN\n");
}
@ -299,21 +278,19 @@ void NEB_spin::run()
if (ulogfile) {
if (verbose) {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
fprintf(ulogfile,"Step MaxReplicaTorque MaxAtomTorque "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
"... ReplicaTorqueN MaxAtomTorqueN\n");
} else {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
fprintf(ulogfile,"Step MaxReplicaTorque MaxAtomTorque "
"GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... "
"RDN PEN\n");
}
}
}
//printf("test before run 4 \n");
print_status();
//printf("test before run 5 \n");
// perform regular NEB_spin for n1steps or until replicas converge
// retrieve PE values from fix NEB_spin and print every nevery iterations
@ -323,35 +300,11 @@ void NEB_spin::run()
timer->init();
timer->barrier_start();
// test import fix_nve scheme
//printf("test 2 atom: i=%d,%g,%g,%g \n",1,sp[1][0],sp[1][1],sp[1][2]);
//error->all(FLERR,"end neb_spin test");
double dts;
while (update->minimize->niter < n1steps) {
//dts = evaluate_dt();
//advance_spins(dts);
//fneb->
//dts = fneb->evaluate_dt();
//fneb->advance_spins(dts);
// no minimizer for spins
update->minimize->run(nevery);
// no minimizer for spins
//update->minimize->run(nevery);
//
// evaluate dts
// loop on spins, damped advance
//
print_status();
if (update->minimize->stop_condition) break;
}
// test neb end
//error->all(FLERR,"end neb_spin test");
timer->barrier_stop();
@ -393,22 +346,19 @@ void NEB_spin::run()
update->minimize->init();
fneb->rclimber = top;
printf("test print 6.2 \n");
update->minimize->setup();
printf("test print 6.3 \n");
if (me_universe == 0) {
if (uscreen) {
if (verbose) {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
fprintf(uscreen,"Step MaxReplicaTorque MaxAtomTorque "
"GradV0 GradV1 GradVc EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN "
"pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(uscreen,"Step MaxReplicaForce MaxAtomForce "
fprintf(uscreen,"Step MaxReplicaTorque MaxAtomTorque "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
@ -416,14 +366,14 @@ void NEB_spin::run()
}
if (ulogfile) {
if (verbose) {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
fprintf(ulogfile,"Step MaxReplicaTorque MaxAtomTorque "
"GradV0 GradV1 GradVc EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN "
"pathangle1 angletangrad1 anglegrad1 gradV1 "
"ReplicaForce1 MaxAtomForce1 pathangle2 angletangrad2 "
"... ReplicaForceN MaxAtomForceN\n");
} else {
fprintf(ulogfile,"Step MaxReplicaForce MaxAtomForce "
fprintf(ulogfile,"Step MaxReplicaTorque MaxAtomTorque "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
@ -441,10 +391,6 @@ void NEB_spin::run()
timer->barrier_start();
while (update->minimize->niter < n2steps) {
//dts = evaluate_dt();
//advance_spins(dts);
//dts = fneb->evaluate_dt();
//fneb->advance_spins(dts);
update->minimize->run(nevery);
print_status();
if (update->minimize->stop_condition) break;
@ -462,174 +408,6 @@ void NEB_spin::run()
update->beginstep = update->endstep = 0;
}
/* ----------------------------------------------------------------------
geodesic distance calculation (Vincenty's formula)
------------------------------------------------------------------------- */
//double NEB_spin::geodesic_distance2(double spi[3], double spj[3])
//{
// double dist;
// double crossx,crossy,crossz;
// double dotx,doty,dotz;
// double crosslen,dots;
//
// crossx = spi[1]*spj[2]-spi[2]*spj[1];
// crossy = spi[2]*spj[0]-spi[0]*spj[2];
// crossz = spi[0]*spj[1]-spi[1]*spj[0];
// crosslen = sqrt(crossx*crossx + crossy*crossy + crossz*crossz);
// dotx = spi[0]*spj[0];
// doty = spi[1]*spj[1];
// dotz = spi[2]*spj[2];
// dots = dotx+doty+dotz;
//
// dist = atan2(crosslen,dots);
//
// return dist;
//}
/* ----------------------------------------------------------------------
evaluate max timestep
---------------------------------------------------------------------- */
//double NEB_spin::evaluate_dt()
//{
// double dtmax;
// double fmsq;
// double fmaxsqone,fmaxsqloc,fmaxsqall;
// int nlocal = atom->nlocal;
// int *mask = atom->mask;
// double **fm = atom->fm;
//
// // finding max fm on this proc.
//
// fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
// for (int i = 0; i < nlocal; i++)
// if (mask[i] & groupbit) {
// fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2];
// fmaxsqone = MAX(fmaxsqone,fmsq);
// }
//
// // finding max fm on this replica
//
// fmaxsqloc = fmaxsqone;
// MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);
//
// // finding max fm over all replicas, if necessary
// // this communicator would be invalid for multiprocess replicas
//
// if (update->multireplica == 1) {
// fmaxsqall = fmaxsqloc;
// MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);
// }
//
// if (fmaxsqall < fmaxsqloc)
// error->all(FLERR,"Incorrect fmaxall calc.");
//
// // define max timestep
// // dividing by 10 the inverse of max frequency
//
// dtmax = MY_2PI/(10.0*sqrt(fmaxsqall));
//
// return dtmax;
//}
/* ----------------------------------------------------------------------
geometric damped advance os spins
---------------------------------------------------------------------- */
//void NEB_spin::advance_spins(double dts)
//{
// //int j=0;
// //int *sametag = atom->sametag;
// int nlocal = atom->nlocal;
// int *mask = atom->mask;
// double **sp = atom->sp;
// double **fm = atom->fm;
// double tdampx,tdampy,tdampz;
// double msq,scale,fm2,energy,dts2;
// double alpha;
// double spi[3],fmi[3];
// double cp[3],g[3];
//
// //cp[0] = cp[1] = cp[2] = 0.0;
// //g[0] = g[1] = g[2] = 0.0;
// dts2 = dts*dts;
//
// // fictitious Gilbert damping of 1
// alpha = 1.0;
//
// // loop on all spins on proc.
//
// if (ireplica != nreplica-1 && ireplica != 0)
// for (int i = 0; i < nlocal; i++)
// if (mask[i] & groupbit) {
//
// spi[0] = sp[i][0];
// spi[1] = sp[i][1];
// spi[2] = sp[i][2];
//
// fmi[0] = fm[i][0];
// fmi[1] = fm[i][1];
// fmi[2] = fm[i][2];
//
// // calc. damping torque
//
// tdampx = -alpha*(fmi[1]*spi[2] - fmi[2]*spi[1]);
// tdampy = -alpha*(fmi[2]*spi[0] - fmi[0]*spi[2]);
// tdampz = -alpha*(fmi[0]*spi[1] - fmi[1]*spi[0]);
//
// // apply advance algorithm (geometric, norm preserving)
//
// fm2 = (tdampx*tdampx+tdampy*tdampy+tdampz*tdampz);
// energy = (sp[i][0]*tdampx)+(sp[i][1]*tdampy)+(sp[i][2]*tdampz);
//
// cp[0] = tdampy*sp[i][2]-tdampz*sp[i][1];
// cp[1] = tdampz*sp[i][0]-tdampx*sp[i][2];
// cp[2] = tdampx*sp[i][1]-tdampy*sp[i][0];
//
// g[0] = sp[i][0]+cp[0]*dts;
// g[1] = sp[i][1]+cp[1]*dts;
// g[2] = sp[i][2]+cp[2]*dts;
//
// g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts2;
// g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts2;
// g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts2;
//
// g[0] /= (1+0.25*fm2*dts2);
// g[1] /= (1+0.25*fm2*dts2);
// g[2] /= (1+0.25*fm2*dts2);
//
// sp[i][0] = g[0];
// sp[i][1] = g[1];
// sp[i][2] = g[2];
//
// // renormalization (check if necessary)
//
// msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
// scale = 1.0/sqrt(msq);
// sp[i][0] *= scale;
// sp[i][1] *= scale;
// sp[i][2] *= scale;
//
// // comm. sp[i] to atoms with same tag (for serial algo)
//
// // no need for simplecticity
// //if (sector_flag == 0) {
// // if (sametag[i] >= 0) {
// // j = sametag[i];
// // while (j >= 0) {
// // sp[j][0] = sp[i][0];
// // sp[j][1] = sp[i][1];
// // sp[j][2] = sp[i][2];
// // j = sametag[j];
// // }
// // }
// //}
// //
//
// }
//}
/* ----------------------------------------------------------------------
read initial config atom coords from file
flag = 0
@ -654,7 +432,6 @@ void NEB_spin::readfile(char *file, int flag)
double xx,yy,zz,delx,dely,delz;
// spin quantities
double musp,spx,spy,spz;
//,delx,dely,delz;
if (me_universe == 0 && screen)
fprintf(screen,"Reading NEB_spin coordinate file(s) ...\n");
@ -703,12 +480,6 @@ void NEB_spin::readfile(char *file, int flag)
double spinit[3],spfinal[3];
int nlocal = atom->nlocal;
// test 1.2
//double **sp = atom->sp;
//printf("test 1.2 atom: i=%d,%g,%g,%g \n",1,sp[1][0],sp[1][1],sp[1][2]);
//error->all(FLERR,"end neb_spin test");
// loop over chunks of lines read from file
// two versions of read_lines_from_file() for world vs universe bcast
// count # of atom coords changed so can check for invalid atom IDs in file
@ -765,44 +536,27 @@ void NEB_spin::readfile(char *file, int flag)
spx = atof(values[5]);
spy = atof(values[6]);
spz = atof(values[7]);
//xx = atof(values[1]);
//yy = atof(values[2]);
//zz = atof(values[3]);
if (flag == 0) {
// here, function interp. spin states
//spinit[0] = x[m][0];
//spinit[1] = x[m][1];
//spinit[2] = x[m][2];
spinit[0] = sp[m][0];
spinit[1] = sp[m][1];
spinit[2] = sp[m][2];
spfinal[0] = spx;
spfinal[1] = spy;
spfinal[2] = spz;
//domain->minimum_image(delx,dely,delz);
// to be used it atomic displacements with pbc
//domain->minimum_image(delx,dely,delz);
// test
//printf("spinit: %g %g %g \n",spinit[0],spinit[1],spinit[2]);
//printf("spfinal bef: %g %g %g \n",spfinal[0],spfinal[1],spfinal[2]);
// interpolate intermediate spin states
initial_rotation(spinit,spfinal,fraction);
// test
//printf("spfinal aft: %g %g %g \n",spfinal[0],spfinal[1],spfinal[2]);
sp[m][0] = spfinal[0];
sp[m][1] = spfinal[1];
sp[m][2] = spfinal[2];
sp[m][3] = musp;
//delx = xx - x[m][0];
//dely = yy - x[m][1];
//delz = zz - x[m][2];
//x[m][0] += fraction*delx;
//x[m][1] += fraction*dely;
//x[m][2] += fraction*delz;
} else {
sp[m][3] = musp;
x[m][0] = xx;
@ -820,12 +574,6 @@ void NEB_spin::readfile(char *file, int flag)
nread += nchunk;
}
// test 1.3
//double **sp = atom->sp;
//printf("test 1.3 atom: i=%d,%g,%g,%g \n",1,sp[1][0],sp[1][1],sp[1][2]);
//error->all(FLERR,"end neb_spin test");
// check that all atom IDs in file were found by a proc
if (flag == 0) {
@ -867,13 +615,11 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
{
// implementing initial rotation using atan2
//atan2(crosslen,dots);
// this is not a sufficient routine,
// we need more accurate verifications
// initial and final and inter ang. values
// initial, final and inter ang. values
double itheta,iphi,ftheta,fphi,ktheta,kphi;
double spix,spiy,spiz,spfx,spfy,spfz;
double spkx,spky,spkz,iknorm;
@ -909,58 +655,6 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
sploc[1] = spky;
sploc[2] = spkz;
//double theta,spdot;
//double inormdot,ispinorm;
//double kix,kiy,kiz;
//double kinorm, ikinorm;
//double crossx,crossy,crossz;
////printf("inside rot, spi %g, spf %g \n",spi[0],sploc[0]);
//spdot = spi[0]*sploc[0]+spi[1]*sploc[1]+spi[2]*sploc[2];
//theta = fraction*acos(spdot);
//printf("inside rot, theta %g \n",theta);
//kix = spi[1]*sploc[2]-spi[2]*sploc[1];
//kiy = spi[2]*sploc[0]-spi[0]*sploc[2];
//kiz = spi[0]*sploc[1]-spi[1]*sploc[0];
//
////printf("inside rot1.1, ki %g %g %g \n",kix,kiy,kiz);
//inormdot = 1.0/sqrt(spdot);
//kinorm = kix*kix+kiy*kiy+kiz*kiz;
//if (kinorm == 0.0) {
// kix = 0.0;
// kiy = 0.0;
// kiz = 0.0;
//} else {
// ikinorm = 1.0/kinorm;
// kix *= ikinorm;
// kiy *= ikinorm;
// kiz *= ikinorm;
//}
////printf("inside rot1.2, kin %g %g %g \n",kix,kiy,kiz);
//crossx = kiy*spi[2]-kiz*spi[1];
//crossy = kiz*spi[0]-kix*spi[2];
//crossz = kix*spi[1]-kiy*spi[0];
//
////printf("inside rot1.3, cross %g %g %g \n",crossx,crossy,crossz);
//sploc[0] = spi[0]*cos(theta)+crossx*sin(theta);
//sploc[1] = spi[1]*cos(theta)+crossy*sin(theta);
//sploc[2] = spi[2]*cos(theta)+crossz*sin(theta);
//
////printf("inside rot2, spf %g %g %g \n",sploc[0],sploc[1],sploc[2]);
//ispinorm = 1.0/sqrt(sploc[0]*sploc[0]+sploc[1]*sploc[1]+sploc[2]*sploc[2]);
//sploc[0] *= ispinorm;
//sploc[1] *= ispinorm;
//sploc[2] *= ispinorm;
//printf("inside rot2, spf %g %g %g \n",sploc[0],sploc[1],sploc[2]);
}
/* ----------------------------------------------------------------------
@ -1007,30 +701,43 @@ void NEB_spin::print_status()
//double fnorm2 = sqrt(update->minimize->fnorm_sqr());
// test fmax spin
int nlocal = atom->nlocal;
double tx,ty,tz;
double tnorm2,local_norm_inf,temp_inf;
double **sp = atom->sp;
double **fm = atom->fm;
double fnorm2;
for (int i = 0; i < nlocal; i++)
fnorm2 += (fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2]);
//for (int i = 0; i < nlocal; i++)
// if (mask[i] & groupbit) {
// fnorm2 += (fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2]);
// }
// calc. magnetic torques
tnorm2 = local_norm_inf = temp_inf = 0.0;
for (int i = 0; i < nlocal; i++) {
tx = (fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]);
ty = (fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]);
tz = (fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]);
tnorm2 += tx*tx + ty*ty + tz*tz;
temp_inf = MAX(fabs(tx),fabs(ty));
temp_inf = MAX(fabs(tz),temp_inf);
local_norm_inf = MAX(temp_inf,local_norm_inf);
}
double fmaxreplica;
MPI_Allreduce(&fnorm2,&fmaxreplica,1,MPI_DOUBLE,MPI_MAX,roots);
//MPI_Allreduce(&fnorm2,&fmaxreplica,1,MPI_DOUBLE,MPI_MAX,roots);
MPI_Allreduce(&tnorm2,&fmaxreplica,1,MPI_DOUBLE,MPI_MAX,roots);
// no minimize->fnorm_inf for spins
//double fnorminf = update->minimize->fnorm_inf();
//double fmaxatom;
//MPI_Allreduce(&fnorminf,&fmaxatom,1,MPI_DOUBLE,MPI_MAX,roots);
double fnorminf = 0.0;
double fmaxatom = 0.0;
MPI_Allreduce(&local_norm_inf,&fnorminf,1,MPI_DOUBLE,MPI_MAX,world);
double fmaxatom;
MPI_Allreduce(&fnorminf,&fmaxatom,1,MPI_DOUBLE,MPI_MAX,roots);
if (verbose) {
freplica = new double[nreplica];
MPI_Allgather(&fnorm2,1,MPI_DOUBLE,&freplica[0],1,MPI_DOUBLE,roots);
//MPI_Allgather(&fnorm2,1,MPI_DOUBLE,&freplica[0],1,MPI_DOUBLE,roots);
MPI_Allgather(&tnorm2,1,MPI_DOUBLE,&freplica[0],1,MPI_DOUBLE,roots);
fmaxatomInRepl = new double[nreplica];
MPI_Allgather(&fnorminf,1,MPI_DOUBLE,&fmaxatomInRepl[0],1,MPI_DOUBLE,roots);
}

6
src/REPLICA/neb_spin.h
View File

@ -45,7 +45,8 @@ class NEB_spin : protected Pointers {
FILE *fp;
int compressed;
double etol; // energy tolerance convergence criterion
double ftol; // force tolerance convergence criterion
//double ftol; // force tolerance convergence criterion
double ttol; // torque tolerance convergence criterion
int n1steps, n2steps; // number of steps in stage 1 and 2
int nevery; // output interval
char *infile; // name of file containing final state
@ -57,9 +58,6 @@ class NEB_spin : protected Pointers {
double *freplica; // force on an image
double *fmaxatomInRepl; // force on an image
//double geodesic_distance2(double *, double *);
//double evaluate_dt();
//void advance_spins(double);
void readfile(char *, int);
void initial_rotation(double *, double *, double);
void open(char *);