Merge pull request #445 from Pakketeretet2/extract_for_bond

Extract for bond
This commit is contained in:
sjplimp 2017-04-11 08:27:31 -06:00 committed by GitHub
commit 06cebb9fb4
7 changed files with 168 additions and 6 deletions

View File

@ -22,6 +22,11 @@ attribute = {pair} or {kspace} or {atom} :l
pparam = parameter to adapt over time
I,J = type pair(s) to set parameter for
v_name = variable with name that calculates value of pparam
{bond} args = bstyle bparam I v_name
bstyle = bond style name, e.g. harmonic
bparam = parameter to adapt over time
I = type bond to set parameter for
v_name = variable with name that calculates value of bparam
{kspace} arg = v_name
v_name = variable with name that calculates scale factor on K-space terms
{atom} args = aparam v_name
@ -42,7 +47,10 @@ keyword = {scale} or {reset} :l
fix 1 all adapt 1 pair soft a 1 1 v_prefactor
fix 1 all adapt 1 pair soft a 2* 3 v_prefactor
fix 1 all adapt 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes
fix 1 all adapt 10 atom diameter v_size :pre
fix 1 all adapt 10 atom diameter v_size
variable ramp_up equal "ramp(0.01,0.5)"
fix stretch all adapt 1 bond harmonic r0 1 v_ramp_up :pre
[Description:]
@ -192,6 +200,19 @@ fix 1 all adapt 1 pair soft a * * v_prefactor :pre
:line
The {bond} keyword uses the specified variable to change the value of
a bond coefficient over time, very similar to how the {pair} keyword
operates. The only difference is that now a bond coefficient for a
given bond type is adapted.
Currently {bond} does not support bond_style hybrid nor bond_style
hybrid/overlay as bond styles. The only bonds that currently are
working with fix_adapt are
"harmonic"_bond_harmonic.html: k,r0: type bonds :tb(c=3,s=:)
:line
The {kspace} keyword used the specified variable as a scale factor on
the energy, forces, virial calculated by whatever K-Space solver is
defined by the "kspace_style"_kspace_style.html command. If the

View File

@ -13,6 +13,7 @@
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "bond_harmonic.h"
#include "atom.h"
#include "neighbor.h"
@ -26,7 +27,10 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondHarmonic::BondHarmonic(LAMMPS *lmp) : Bond(lmp) {}
BondHarmonic::BondHarmonic(LAMMPS *lmp) : Bond(lmp)
{
reinitflag = 1;
}
/* ---------------------------------------------------------------------- */
@ -196,3 +200,16 @@ double BondHarmonic::single(int type, double rsq, int i, int j,
if (r > 0.0) fforce = -2.0*rk/r;
return rk*dr;
}
/* ----------------------------------------------------------------------
Return ptr to internal members upon request.
------------------------------------------------------------------------ */
void *BondHarmonic::extract( char *str, int &dim )
{
dim = 1;
if( strcmp(str,"kappa")==0) return (void*) k;
if( strcmp(str,"r0")==0) return (void*) r0;
return NULL;
}

View File

@ -36,6 +36,7 @@ class BondHarmonic : public Bond {
void read_restart(FILE *);
void write_data(FILE *);
double single(int, double, int, int, double &);
virtual void *extract(char *, int &);
protected:
double *k,*r0;

View File

@ -292,3 +292,14 @@ double Bond::memory_usage()
bytes += comm->nthreads*maxvatom*6 * sizeof(double);
return bytes;
}
/* -----------------------------------------------------------------------
Reset all type-based bond params via init.
-------------------------------------------------------------------------- */
void Bond::reinit()
{
if (!reinitflag)
error->all(FLERR,"Fix adapt interface to this bond style not supported");
init();
}

View File

@ -30,6 +30,8 @@ class Bond : protected Pointers {
double virial[6]; // accumulated virial
double *eatom,**vatom; // accumulated per-atom energy/virial
int reinitflag; // 1 if compatible with fix adapt and alike
// KOKKOS host/device flag and data masks
ExecutionSpace execution_space;
@ -49,6 +51,8 @@ class Bond : protected Pointers {
virtual void write_data(FILE *) {}
virtual double single(int, double, int, int, double &) = 0;
virtual double memory_usage();
virtual void *extract(char *, int &) {return NULL;}
virtual void reinit();
void write_file(int, char**);

View File

@ -16,6 +16,7 @@
#include <stdlib.h>
#include "fix_adapt.h"
#include "atom.h"
#include "bond.h"
#include "update.h"
#include "group.h"
#include "modify.h"
@ -35,7 +36,7 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{PAIR,KSPACE,ATOM};
enum{PAIR,KSPACE,ATOM,BOND};
enum{DIAMETER,CHARGE};
/* ---------------------------------------------------------------------- */
@ -68,6 +69,10 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL)
if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command");
nadapt++;
iarg += 3;
} else if (strcmp(arg[iarg],"bond") == 0 ){
if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command");
nadapt++;
iarg += 5;
} else break;
}
@ -103,6 +108,25 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL)
} else error->all(FLERR,"Illegal fix adapt command");
nadapt++;
iarg += 6;
} else if (strcmp(arg[iarg],"bond") == 0 ){
if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command");
adapt[nadapt].which = BOND;
int n = strlen(arg[iarg+1]) + 1;
adapt[nadapt].bstyle = new char[n];
strcpy(adapt[nadapt].bstyle,arg[iarg+1]);
n = strlen(arg[iarg+2]) + 1;
adapt[nadapt].bparam = new char[n];
adapt[nadapt].bond = NULL;
strcpy(adapt[nadapt].bparam,arg[iarg+2]);
force->bounds(FLERR,arg[iarg+3],atom->ntypes,
adapt[nadapt].ilo,adapt[nadapt].ihi);
if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) {
n = strlen(&arg[iarg+4][2]) + 1;
adapt[nadapt].var = new char[n];
strcpy(adapt[nadapt].var,&arg[iarg+4][2]);
} else error->all(FLERR,"Illegal fix adapt command");
nadapt++;
iarg += 5;
} else if (strcmp(arg[iarg],"kspace") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command");
adapt[nadapt].which = KSPACE;
@ -160,6 +184,13 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL)
for (int m = 0; m < nadapt; m++)
if (adapt[m].which == PAIR)
memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
// allocate bond style arrays:
n = atom->nbondtypes;
for (int m = 0; m < nadapt; ++m)
if (adapt[m].which == BOND)
// For now just use same storage and fake it to be one-dimensional:
memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
}
/* ---------------------------------------------------------------------- */
@ -172,6 +203,10 @@ FixAdapt::~FixAdapt()
delete [] adapt[m].pstyle;
delete [] adapt[m].pparam;
memory->destroy(adapt[m].array_orig);
} else if (adapt[m].which == BOND) {
delete [] adapt[m].bstyle;
delete [] adapt[m].bparam;
memory->destroy(adapt[m].vector_orig);
}
}
delete [] adapt;
@ -282,6 +317,7 @@ void FixAdapt::init()
// setup and error checks
anypair = 0;
anybond = 0;
for (int m = 0; m < nadapt; m++) {
Adapt *ad = &adapt[m];
@ -350,7 +386,44 @@ void FixAdapt::init()
}
delete [] pstyle;
} else if (ad->which == BOND){
ad->bond = NULL;
anybond = 1;
// Use same routines from pair to strip any suffices:
int n = strlen(ad->bstyle) + 1;
char *bstyle = new char[n];
strcpy(bstyle,ad->bstyle);
if (lmp->suffix_enable) {
int len = 2 + strlen(bstyle) + strlen(lmp->suffix);
char *bsuffix = new char[len];
strcpy(bsuffix,bstyle);
strcat(bsuffix,"/");
strcat(bsuffix,lmp->suffix);
ad->bond = force->bond_match(bsuffix);
delete [] bsuffix;
}
// If not set grab regular one instead:
if (ad->bond == NULL) ad->bond = force->bond_match(bstyle);
if (ad->bond == NULL )
error->all(FLERR,"Fix adapt bond style does not exist");
void *ptr = ad->bond->extract(ad->bparam,ad->bdim);
if (ptr == NULL)
error->all(FLERR,"Fix adapt bond style param not supported");
// For bond styles you should use a vector
if (ad->bdim == 1) ad->vector = (double *) ptr;
if (strcmp(force->bond_style,"hybrid") == 0 ||
strcmp(force->bond_style,"hybrid_overlay") == 0)
error->all(FLERR,"Fix adapt does not support bond_style hybrid");
delete [] bstyle;
} else if (ad->which == KSPACE) {
if (force->kspace == NULL)
error->all(FLERR,"Fix adapt kspace style does not exist");
@ -368,7 +441,7 @@ void FixAdapt::init()
}
}
// make copy of original pair array values
// make copy of original pair/bond array values
for (int m = 0; m < nadapt; m++) {
Adapt *ad = &adapt[m];
@ -378,7 +451,12 @@ void FixAdapt::init()
ad->array_orig[i][j] = ad->array[i][j];
}else if (ad->which == PAIR && ad->pdim == 0){
ad->scalar_orig = *ad->scalar;
}else if (ad->which == BOND && ad->bdim == 1){
for (i = ad->ilo; i <= ad->ihi; ++i )
ad->vector_orig[i] = ad->vector[i];
}
}
// fixes that store initial per-atom values
@ -470,6 +548,18 @@ void FixAdapt::change_settings()
ad->array[i][j] = value;
}
// set bond type array values:
} else if (ad->which == BOND) {
if (ad->bdim == 1){
if (scaleflag)
for (i = ad->ilo; i <= ad->ihi; ++i )
ad->vector[i] = value*ad->vector_orig[i];
else
for (i = ad->ilo; i <= ad->ihi; ++i )
ad->vector[i] = value;
}
// set kspace scale factor
} else if (ad->which == KSPACE) {
@ -532,6 +622,14 @@ void FixAdapt::change_settings()
}
}
}
if (anybond) {
for (int m = 0; m < nadapt; ++m ) {
Adapt *ad = &adapt[m];
if (ad->which == BOND) {
ad->bond->reinit();
}
}
}
// reset KSpace charges if charges have changed
@ -554,6 +652,12 @@ void FixAdapt::restore_settings()
ad->array[i][j] = ad->array_orig[i][j];
}
} else if (ad->which == BOND) {
if (ad->pdim == 1) {
for (int i = ad->ilo; i <= ad->ihi; i++)
ad->vector[i] = ad->vector_orig[i];
}
} else if (ad->which == KSPACE) {
*kspace_scale = 1.0;
@ -588,6 +692,7 @@ void FixAdapt::restore_settings()
}
if (anypair) force->pair->reinit();
if (anybond) force->bond->reinit();
if (chgflag && force->kspace) force->kspace->qsum_qsq();
}

View File

@ -43,7 +43,7 @@ class FixAdapt : public Fix {
private:
int nadapt,resetflag,scaleflag;
int anypair;
int anypair, anybond;
int nlevels_respa;
char *id_fix_diam,*id_fix_chg;
class FixStore *fix_diam,*fix_chg;
@ -52,12 +52,15 @@ class FixAdapt : public Fix {
int which,ivar;
char *var;
char *pstyle,*pparam;
char *bstyle,*bparam;
int ilo,ihi,jlo,jhi;
int pdim;
int pdim,bdim;
double *scalar,scalar_orig;
double *vector,*vector_orig;
double **array,**array_orig;
int aparam;
class Pair *pair;
class Bond *bond;
};
Adapt *adapt;