/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, [email protected]
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory. **/
#include "fix_area.h"
#include <mpi.h>
#include <cstring>
#include <cstdlib>
#include "atom.h"
#include "atom_masks.h"
#include "comm.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
#include "force.h"
#include "atom_vec.h"
#include "molecule.h"
using namespace LAMMPS_NS;
using namespace FixConst;
FixArea::FixArea(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),atml(nullptr),area(nullptr),potarea(nullptr)
{
if (strcmp(style,"area") != 0 && narg < 5)
error->all(FLERR,"Illegal fix dumbbell command: not enough args");
A0 = utils::numeric(FLERR,arg[3],false,lmp);
Xi =utils::numeric(FLERR,arg[4],false,lmp);
n =utils::numeric(FLERR,arg[5],false,lmp); //number of atoms per molecule
memory->create(atml, n_mol +1, n , "FixArea:atml");
memory->create(area,n_mol+1,"FixArea:area");
memory->create(potarea,n_mol+1,"FixArea:potarea");
}
int FixArea::setmask()
{
int mask = 0;
mask |= POST_FORCE;
return mask;
}
void FixArea::post_force(int /*vflag*/)
{
double fx, fy ;
tagint *molecule = atom->molecule;
int id_mol = 0 ;
double **x = atom->x;
double **f = atom->f;
int *molatom = atom->molatom;
int molecular = atom->molecular;
imageint *image = atom->image;
int jatom;
for (int i = 0; i < atom->nlocal; i++) {
// molecule id starts from 1. max(id_mol) equals to the number of molecules in the system
id_mol = std::max(id_mol, molecule[i]);
}
MPI_Allreduce(&id_mol, &n_mol, 1, MPI_LMP_TAGINT, MPI_MAX, world);
int **atml_tmp;
memory->create(atml_tmp, n_mol +1, n , "FixArea:atml_tmp");
memset(&atml_tmp[0][0], 0, sizeof(int) * n_mol +1 * n+1);
if(molecular == 1){
for (int i = 1; i < n_mol+1; ++i) {
for (int k = 0; k < atom->nlocal; ++k) {
if (molecule[k] == i) {
jatom = molatom[k];
atml_tmp[i][jatom] = k;
}
}
}
}
else{
error->all(FLERR,"not a molecular template system ");
}
memory->destroy(atml);
memory->create(atml, n_mol +1, n +1, "FixArea:atml");
MPI_Allreduce(atml_tmp,atml, (n_mol + 1)*(n+1), MPI_INT, MPI_SUM, world);
memory->destroy(atml_tmp);
memory->destroy(area);
memory->create(area,n_mol+1,"FixArea:area");
for(int i = 1 ; i< n_mol + 1;++i){
double unwrap1[3];
double unwrap2[3];
double unwrap3[3];
double unwrap4[3];
for(int jj=1 ; jj< n+1; ++jj){
domain->unmap(x[atml[i][jj]], image[atml[i][jj]], unwrap1);
domain->unmap(x[atml[i][jj+1]], image[atml[i][jj+1]], unwrap2);
area[i] += unwrap1[0]*unwrap2[1] - unwrap2[0] * unwrap1[0];
}
domain->unmap(x[atml[i][n]], image[atml[i][n]], unwrap3);
domain->unmap(x[atml[i][1]], image[atml[i][1]], unwrap4);
area[i] += unwrap3[0]*unwrap4[1]-unwrap3[1] * unwrap4[0];
fprintf(screen, "area = %f \n", area[i]);
}
memory->destroy(potarea);
memory->create(potarea,n_mol+1,"FixArea:potarea");
for (int i = 1 ; i< n_mol+1; ++i){
potarea[i]= 0.5 * Xi *(1-(area[i]/A0))*(1-(area[i]/A0));
}
for(int i=1 ;i < n_mol+1; i++){
for(int j=1 ; j< n+1;++j){
if(j == 1){
fx= Xi *A0 *(1-(area[molecule[i]]/A0))*(x[atml[i][j+1]][1]-x[atml[i][n]][1]);
fy= Xi *A0 *(1-(area[molecule[i]]/A0))*(x[atml[i][n]][0]-x[atml[i][j+1]][0]);
f[atml[i][j]][0] += fx; // unit vector rotated CW
f[atml[i][j]][1] += fy ;
}
if(j == n){
fx= -Xi *A0 *(1-(area[molecule[i]]/A0))*(x[atml[i][1]][1]-x[atml[i][j-1]][1]);
fy= -Xi *A0 *(1-(area[molecule[i]]/A0))*(x[atml[i][j-1]][0]-x[atml[i][1]][0]);
f[atml[i][j]][0] += fx;
f[atml[i][j]][1] += fy ;
}
else{
fx= -Xi *A0 *(1-(area[molecule[i]]/A0))*(x[atml[i][j+1]][1]-x[atml[i][j-1]][1]);
fy= -Xi *A0 *(1-(area[molecule[i]]/A0))*(x[atml[i][j-1]][0]-x[atml[i][j +1]][0]);
f[atml[i][j]][0] += fx;
f[atml[i][j]][1] += fy ;
}
}
}
}
I am writing a custom fix whih calculate the potential energy of a ring molecule associated with its area ,and the corresponding force due to it ,
while compilation i get the above error .The above is my fix_area.cpp file