Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85646639
lmpqst.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Sep 30, 16:43
Size
7 KB
Mime Type
text/x-c
Expires
Wed, Oct 2, 16:43 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21075855
Attached To
rLAMMPS lammps
lmpqst.cpp
View Options
// lmpqst = umbrella driver to couple LAMMPS + Quest
// for MD using quantum forces
// Syntax: lmpqst Niter in.lammps in.quest
// Niter = # of MD iterations
// in.lammps = LAMMPS input script
// in.quest = Quest input script
#include "mpi.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "stdint.h"
#include "many2one.h"
#include "one2many.h"
#include "files.h"
#include "memory.h"
#include "error.h"
#define QUOTE_(x) #x
#define QUOTE(x) QUOTE_(x)
#include "lmppath.h"
#include QUOTE(LMPPATH/src/lammps.h)
#include QUOTE(LMPPATH/src/library.h)
#include QUOTE(LMPPATH/src/input.h)
#include QUOTE(LMPPATH/src/modify.h)
#include QUOTE(LMPPATH/src/fix.h)
#include QUOTE(LMPPATH/src/fix_external.h)
#include "qstexe.h"
using namespace LAMMPS_NS;
#define ANGSTROM_per_BOHR 0.529
#define EV_per_RYDBERG 13.6056923
void quest_callback(void *, bigint, int, int *, double **, double **);
struct Info {
int me;
Memory *memory;
LAMMPS *lmp;
char *quest_input;
};
/* ---------------------------------------------------------------------- */
int main(int narg, char **arg)
{
int n;
char str[128];
// setup MPI
MPI_Init(&narg,&arg);
MPI_Comm comm = MPI_COMM_WORLD;
int me,nprocs;
MPI_Comm_rank(comm,&me);
MPI_Comm_size(comm,&nprocs);
Memory *memory = new Memory(comm);
Error *error = new Error(comm);
// command-line args
if (narg != 4) error->all("Syntax: lmpqst Niter in.lammps in.quest");
int niter = atoi(arg[1]);
n = strlen(arg[2]) + 1;
char *lammps_input = new char[n];
strcpy(lammps_input,arg[2]);
n = strlen(arg[3]) + 1;
char *quest_input = new char[n];
strcpy(quest_input,arg[3]);
// instantiate LAMMPS
LAMMPS *lmp = new LAMMPS(0,NULL,MPI_COMM_WORLD);
// create simulation in LAMMPS from in.lammps
lmp->input->file(lammps_input);
// make info avaiable to callback function
Info info;
info.me = me;
info.memory = memory;
info.lmp = lmp;
info.quest_input = quest_input;
// set callback to Quest inside fix external
// this could also be done thru Python, using a ctypes callback
int ifix = lmp->modify->find_fix("2");
FixExternal *fix = (FixExternal *) lmp->modify->fix[ifix];
fix->set_callback(quest_callback,&info);
// run LAMMPS for Niter
// each time it needs forces, it will invoke quest_callback
sprintf(str,"run %d",niter);
lmp->input->one(str);
// clean up
delete lmp;
delete memory;
delete error;
delete [] lammps_input;
delete [] quest_input;
MPI_Finalize();
}
/* ----------------------------------------------------------------------
callback to Quest with atom IDs and coords from each proc
invoke Quest to compute forces, load them into f for LAMMPS to use
f can be NULL if proc owns no atoms
------------------------------------------------------------------------- */
void quest_callback(void *ptr, bigint ntimestep,
int nlocal, int *id, double **x, double **f)
{
int i,j;
char str[128];
Info *info = (Info *) ptr;
// boxlines = LAMMPS box size converted into Quest lattice vectors
char **boxlines = NULL;
if (info->me == 0) {
boxlines = new char*[3];
for (i = 0; i < 3; i++) boxlines[i] = new char[128];
}
double boxxlo = *((double *) lammps_extract_global(info->lmp,"boxxlo"));
double boxxhi = *((double *) lammps_extract_global(info->lmp,"boxxhi"));
double boxylo = *((double *) lammps_extract_global(info->lmp,"boxylo"));
double boxyhi = *((double *) lammps_extract_global(info->lmp,"boxyhi"));
double boxzlo = *((double *) lammps_extract_global(info->lmp,"boxzlo"));
double boxzhi = *((double *) lammps_extract_global(info->lmp,"boxzhi"));
double boxxy = *((double *) lammps_extract_global(info->lmp,"xy"));
double boxxz = *((double *) lammps_extract_global(info->lmp,"xz"));
double boxyz = *((double *) lammps_extract_global(info->lmp,"yz"));
double xprd = (boxxhi-boxxlo)/ANGSTROM_per_BOHR;
double yprd = (boxyhi-boxylo)/ANGSTROM_per_BOHR;
double zprd = (boxzhi-boxzlo)/ANGSTROM_per_BOHR;
double xy = boxxy/ANGSTROM_per_BOHR;
double xz = boxxz/ANGSTROM_per_BOHR;
double yz = boxyz/ANGSTROM_per_BOHR;
if (info->me == 0) {
sprintf(boxlines[0],"%g %g %g\n",xprd,0.0,0.0);
sprintf(boxlines[1],"%g %g %g\n",xy,yprd,0.0);
sprintf(boxlines[2],"%g %g %g\n",xz,yz,zprd);
}
// xlines = x for atoms on each proc converted to text lines
// xlines is suitable for insertion into Quest input file
// convert LAMMPS Angstroms to Quest bohr
int natoms;
MPI_Allreduce(&nlocal,&natoms,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
Many2One *lmp2qst = new Many2One(MPI_COMM_WORLD);
lmp2qst->setup(nlocal,id,natoms);
char **xlines = NULL;
double **xquest = NULL;
if (info->me == 0) {
xquest = info->memory->create_2d_double_array(natoms,3,"lmpqst:xquest");
xlines = new char*[natoms];
for (i = 0; i < natoms; i++) xlines[i] = new char[128];
}
if (info->me == 0) lmp2qst->gather(&x[0][0],3,&xquest[0][0]);
else lmp2qst->gather(&x[0][0],3,NULL);
if (info->me == 0) {
for (i = 0; i < natoms; i++) {
xquest[i][0] /= ANGSTROM_per_BOHR;
xquest[i][1] /= ANGSTROM_per_BOHR;
xquest[i][2] /= ANGSTROM_per_BOHR;
}
for (i = 0; i < natoms; i++) {
sprintf(xlines[i],"%d %d %g %g %g\n",i+1,1,
xquest[i][0],xquest[i][1],xquest[i][2]);
}
}
// one-processor tasks:
// whack all lcao.* files
// cp quest_input to lcao.in
// replace atom coords section of lcao.in with new atom coords
// run Quest on one proc, save screen output to file
// flines = atom forces extracted from Quest screen file
// fquest = atom forces
// convert Quest Ryd/bohr to LAMMPS eV/Angstrom
char **flines = NULL;
double **fquest = NULL;
if (info->me == 0) {
fquest = info->memory->create_2d_double_array(natoms,3,"lmpqst:fquest");
flines = new char*[natoms];
for (i = 0; i < natoms; i++) flines[i] = new char[128];
}
if (info->me == 0) {
system("rm lcao.*");
sprintf(str,"cp %s lcao.in",info->quest_input);
system(str);
sprintf(str,"cp %s lcao.x",QUOTE(QUEST));
system(str);
replace("lcao.in","primitive lattice vectors",3,boxlines);
replace("lcao.in","atom, type, position vector",natoms,xlines);
system("lcao.x > lcao.screen");
extract("lcao.screen","atom x force "
"y force z force",natoms,flines);
int itmp;
for (i = 0; i < natoms; i++)
sscanf(flines[i],"%d %lg %lg %lg",&itmp,
&fquest[i][0],&fquest[i][1],&fquest[i][2]);
for (i = 0; i < natoms; i++) {
fquest[i][0] *= EV_per_RYDBERG / ANGSTROM_per_BOHR;
fquest[i][1] *= EV_per_RYDBERG / ANGSTROM_per_BOHR;
fquest[i][2] *= EV_per_RYDBERG / ANGSTROM_per_BOHR;
}
}
// convert fquest on one proc into f for atoms on each proc
One2Many *qst2lmp = new One2Many(MPI_COMM_WORLD);
qst2lmp->setup(natoms,nlocal,id);
double *fvec = NULL;
if (f) fvec = &f[0][0];
if (info->me == 0) qst2lmp->scatter(&fquest[0][0],3,fvec);
else qst2lmp->scatter(NULL,3,fvec);
// clean up
// some data only exists on proc 0
delete lmp2qst;
delete qst2lmp;
info->memory->destroy_2d_double_array(xquest);
info->memory->destroy_2d_double_array(fquest);
if (boxlines) {
for (i = 0; i < 3; i++) delete [] boxlines[i];
delete [] boxlines;
}
if (xlines) {
for (i = 0; i < natoms; i++) delete [] xlines[i];
delete [] xlines;
}
if (flines) {
for (i = 0; i < natoms; i++) delete [] flines[i];
delete [] flines;
}
}
Event Timeline
Log In to Comment