Commit JT 032619

- finish merge of min_spin.cpp
- test output
This commit is contained in:
julient31 2019-03-26 12:51:51 -06:00
parent b8747ce8fb
commit 3cabfd13cd
2 changed files with 2 additions and 117 deletions

View File

@ -51,11 +51,7 @@ MinSpin::MinSpin(LAMMPS *lmp) : Min(lmp) {}
void MinSpin::init()
{
alpha_damp = 1.0;
<<<<<<< HEAD
discret_factor = 10.0;
=======
discrete_factor = 10.0;
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
Min::init();
@ -88,15 +84,9 @@ int MinSpin::modify_param(int narg, char **arg)
alpha_damp = force->numeric(FLERR,arg[1]);
return 2;
}
<<<<<<< HEAD
if (strcmp(arg[0],"discret_factor") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
discret_factor = force->numeric(FLERR,arg[1]);
=======
if (strcmp(arg[0],"discrete_factor") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
discrete_factor = force->numeric(FLERR,arg[1]);
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
return 2;
}
return 0;
@ -114,17 +104,10 @@ void MinSpin::reset_vectors()
// size sp is 4N vector
nvec = 4 * atom->nlocal;
if (nvec) spvec = atom->sp[0];
<<<<<<< HEAD
nvec = 3 * atom->nlocal;
if (nvec) fmvec = atom->fm[0];
=======
nvec = 3 * atom->nlocal;
if (nvec) fmvec = atom->fm[0];
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
if (nvec) xvec = atom->x[0];
if (nvec) fvec = atom->f[0];
}
@ -136,11 +119,7 @@ void MinSpin::reset_vectors()
int MinSpin::iterate(int maxiter)
{
bigint ntimestep;
<<<<<<< HEAD
double fmdotfm,fmdotfmall;
=======
double fmdotfm;
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
int flag,flagall;
for (int iter = 0; iter < maxiter; iter++) {
@ -153,21 +132,12 @@ int MinSpin::iterate(int maxiter)
// optimize timestep accross processes / replicas
// need a force calculation for timestep optimization
<<<<<<< HEAD
energy_force(0);
dts = evaluate_dt();
// apply damped precessional dynamics to the spins
=======
energy_force(0);
dts = evaluate_dt();
// apply damped precessional dynamics to the spins
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
advance_spins(dts);
eprevious = ecurrent;
@ -230,18 +200,10 @@ double MinSpin::evaluate_dt()
double fmsq;
double fmaxsqone,fmaxsqloc,fmaxsqall;
int nlocal = atom->nlocal;
<<<<<<< HEAD
int *mask = atom->mask;
double **fm = atom->fm;
// finding max fm on this proc.
=======
double **fm = atom->fm;
// finding max fm on this proc.
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
for (int i = 0; i < nlocal; i++) {
fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2];
@ -249,17 +211,10 @@ double MinSpin::evaluate_dt()
}
// finding max fm on this replica
<<<<<<< HEAD
fmaxsqloc = fmaxsqone;
MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);
=======
fmaxsqloc = fmaxsqone;
MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
// finding max fm over all replicas, if necessary
// this communicator would be invalid for multiprocess replicas
@ -272,17 +227,10 @@ double MinSpin::evaluate_dt()
if (fmaxsqall == 0.0)
error->all(FLERR,"Incorrect fmaxsqall calculation");
<<<<<<< HEAD
// define max timestep by dividing by the
// inverse of max frequency by discret_factor
dtmax = MY_2PI/(discret_factor*sqrt(fmaxsqall));
=======
// define max timestep by dividing by the
// inverse of max frequency by discrete_factor
dtmax = MY_2PI/(discrete_factor*sqrt(fmaxsqall));
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
return dtmax;
}
@ -294,51 +242,17 @@ double MinSpin::evaluate_dt()
void MinSpin::advance_spins(double dts)
{
int nlocal = atom->nlocal;
<<<<<<< HEAD
int *mask = atom->mask;
=======
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
double **sp = atom->sp;
double **fm = atom->fm;
double tdampx,tdampy,tdampz;
double msq,scale,fm2,energy,dts2;
double cp[3],g[3];
<<<<<<< HEAD
dts2 = dts*dts;
=======
dts2 = dts*dts;
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
// loop on all spins on proc.
for (int i = 0; i < nlocal; i++) {
<<<<<<< HEAD
// calc. damping torque
tdampx = -alpha_damp*(fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]);
tdampy = -alpha_damp*(fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]);
tdampz = -alpha_damp*(fm[i][0]*sp[i][1] - fm[i][1]*sp[i][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] += (tdampx*energy-0.5*sp[i][0]*fm2)*0.5*dts2;
g[1] += (tdampy*energy-0.5*sp[i][1]*fm2)*0.5*dts2;
g[2] += (tdampz*energy-0.5*sp[i][2]*fm2)*0.5*dts2;
=======
// calc. damping torque
@ -363,34 +277,22 @@ void MinSpin::advance_spins(double dts)
g[1] += (tdampy*energy-0.5*sp[i][1]*fm2)*0.5*dts2;
g[2] += (tdampz*energy-0.5*sp[i][2]*fm2)*0.5*dts2;
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
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];
<<<<<<< HEAD
sp[i][2] = g[2];
// renormalization (check if necessary)
=======
sp[i][2] = g[2];
// renormalization (check if necessary)
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
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;
<<<<<<< HEAD
=======
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
// no comm. to atoms with same tag
// because no need for simplecticity
}
@ -402,23 +304,13 @@ void MinSpin::advance_spins(double dts)
double MinSpin::fmnorm_sqr()
{
<<<<<<< HEAD
int i,n;
double *fmatom;
=======
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
int nlocal = atom->nlocal;
double tx,ty,tz;
double **sp = atom->sp;
double **fm = atom->fm;
// calc. magnetic torques
<<<<<<< HEAD
=======
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
double local_norm2_sqr = 0.0;
for (int i = 0; i < nlocal; i++) {
tx = (fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]);
@ -427,19 +319,11 @@ double MinSpin::fmnorm_sqr()
local_norm2_sqr += tx*tx + ty*ty + tz*tz;
}
<<<<<<< HEAD
// no extra atom calc. for spins
if (nextra_atom)
error->all(FLERR,"extra atom option not available yet");
=======
// no extra atom calc. for spins
if (nextra_atom)
error->all(FLERR,"extra atom option not available yet");
>>>>>>> e2e4fe2cf7a95a04ae6a7de93d9b72ad56f0b620
double norm2_sqr = 0.0;
MPI_Allreduce(&local_norm2_sqr,&norm2_sqr,1,MPI_DOUBLE,MPI_SUM,world);

View File

@ -901,7 +901,8 @@ void NEB_spin::print_status()
gradvnorm0,gradvnorm1,gradvnormc);
fprintf(ulogfile,"%12.8g %12.8g %12.8g ",ebf,ebr,endpt);
for (int i = 0; i < nreplica; i++)
fprintf(ulogfile,"%12.8g %12.8g ",rdist[i],all[i][0]);
//fprintf(ulogfile,"%12.8g %12.8g ",rdist[i],all[i][0]);
fprintf(ulogfile,"%12.8g %12.8g %12.8g %12.8g ",rdist[i],all[i][0],all[i][2],all[i][5]);
if (verbose) {
fprintf(ulogfile,"%12.5g %12.5g %12.5g %12.5g %12.5g %12.5g",
NAN,180-acos(all[0][5])*todeg,180-acos(all[0][6])*todeg,