1 /* ----------------------------------------------------------------------
2    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3    https://www.lammps.org/
4    Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
5 
6    Copyright (2003) Sandia Corporation.  Under the terms of Contract
7    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8    certain rights in this software.  This software is distributed under
9    the GNU General Public License.
10 
11    See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13 
14 // c++_driver = simple example of how an umbrella program
15 //              can invoke LAMMPS as a library on some subset of procs
16 // Syntax: simpleCC P in.lammps
17 //         P = # of procs to run LAMMPS on
18 //             must be <= # of procs the driver code itself runs on
19 //         in.lammps = LAMMPS input script
20 // See README for compilation instructions
21 
22 #include <cstdio>
23 #include <cstdlib>
24 #include <cstring>
25 #include <mpi.h>
26 
27 // these are LAMMPS include files
28 #include "lammps.h"
29 #include "input.h"
30 #include "atom.h"
31 #include "library.h"
32 
33 using namespace LAMMPS_NS;
34 
main(int narg,char ** arg)35 int main(int narg, char **arg)
36 {
37   // setup MPI and various communicators
38   // driver runs on all procs in MPI_COMM_WORLD
39   // comm_lammps only has 1st P procs (could be all or any subset)
40 
41   MPI_Init(&narg,&arg);
42 
43   if (narg != 3) {
44     printf("Syntax: simpleCC P in.lammps\n");
45     exit(1);
46   }
47 
48   int me,nprocs;
49   MPI_Comm_rank(MPI_COMM_WORLD,&me);
50   MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
51 
52   int nprocs_lammps = atoi(arg[1]);
53   if (nprocs_lammps > nprocs) {
54     if (me == 0)
55       printf("ERROR: LAMMPS cannot use more procs than available\n");
56     MPI_Abort(MPI_COMM_WORLD,1);
57   }
58 
59   int lammps;
60   if (me < nprocs_lammps) lammps = 1;
61   else lammps = MPI_UNDEFINED;
62   MPI_Comm comm_lammps;
63   MPI_Comm_split(MPI_COMM_WORLD,lammps,0,&comm_lammps);
64 
65   // open LAMMPS input script
66 
67   FILE *fp;
68   if (me == 0) {
69     fp = fopen(arg[2],"r");
70     if (fp == NULL) {
71       printf("ERROR: Could not open LAMMPS input script\n");
72       MPI_Abort(MPI_COMM_WORLD,1);
73     }
74   }
75 
76   // run the input script thru LAMMPS one line at a time until end-of-file
77   // driver proc 0 reads a line, Bcasts it to all procs
78   // (could just send it to proc 0 of comm_lammps and let it Bcast)
79   // all LAMMPS procs call input->one() on the line
80 
81   LAMMPS *lmp = NULL;
82   if (lammps == 1) lmp = new LAMMPS(0,NULL,comm_lammps);
83 
84   int n;
85   char line[1024];
86   while (1) {
87     if (me == 0) {
88       if (fgets(line,1024,fp) == NULL) n = 0;
89       else n = strlen(line) + 1;
90       if (n == 0) fclose(fp);
91     }
92     MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
93     if (n == 0) break;
94     MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD);
95     if (lammps == 1) lammps_command(lmp,line);
96   }
97 
98   // run 10 more steps
99   // get coords from LAMMPS
100   // change coords of 1st atom
101   // put coords back into LAMMPS
102   // run a single step with changed coords
103 
104   double *x = NULL;
105   double *v = NULL;
106 
107   if (lammps == 1) {
108     lmp->input->one("run 10");
109 
110     int natoms = static_cast<int> (lmp->atom->natoms);
111     x = new double[3*natoms];
112     v = new double[3*natoms];
113     lammps_gather_atoms(lmp,(char *) "x",1,3,x);
114     lammps_gather_atoms(lmp,(char *) "v",1,3,v);
115     double epsilon = 0.1;
116     x[0] += epsilon;
117     lammps_scatter_atoms(lmp,(char *) "x",1,3,x);
118 
119     // these 2 lines are the same
120 
121     // lammps_command(lmp,"run 1");
122     lmp->input->one("run 1");
123   }
124 
125   // extract force on single atom two different ways
126 
127   if (lammps == 1) {
128     double **f = (double **) lammps_extract_atom(lmp,(char *) "f");
129     printf("Force on 1 atom via extract_atom: %g\n",f[0][0]);
130 
131     double *fx = (double *)
132       lammps_extract_variable(lmp,(char *) "fx",(char *) "all");
133     printf("Force on 1 atom via extract_variable: %g\n",fx[0]);
134   }
135 
136   // use commands_string() and commands_list() to invoke more commands
137 
138   const char *strtwo = (char *) "run 10\nrun 20";
139   if (lammps == 1) lammps_commands_string(lmp,strtwo);
140 
141   const char *cmds[2];
142   cmds[0] = (char *) "run 10";
143   cmds[1] = (char *) "run 20";
144   if (lammps == 1) lammps_commands_list(lmp,2,cmds);
145 
146   // delete all atoms
147   // create_atoms() to create new ones with old coords, vels
148   // initial thermo should be same as step 20
149 
150   int *type = NULL;
151 
152   if (lammps == 1) {
153     int natoms = static_cast<int> (lmp->atom->natoms);
154     type = new int[natoms];
155     for (int i = 0; i < natoms; i++) type[i] = 1;
156 
157     lmp->input->one("delete_atoms group all");
158     lammps_create_atoms(lmp,natoms,NULL,type,x,v,NULL,0);
159     lmp->input->one("run 10");
160   }
161 
162   delete [] x;
163   delete [] v;
164   delete [] type;
165 
166   // close down LAMMPS
167 
168   delete lmp;
169 
170   // close down MPI
171 
172   if (lammps == 1) MPI_Comm_free(&comm_lammps);
173   MPI_Barrier(MPI_COMM_WORLD);
174   MPI_Finalize();
175 }
176