1 // clang-format off
2 /* -*- c++ -*- ----------------------------------------------------------
3    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4    https://www.lammps.org/, Sandia National Laboratories
5    Steve Plimpton, sjplimp@sandia.gov
6 
7    Copyright (2003) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14 
15 /* ----------------------------------------------------------------------
16    This file is part of the MGPT implementation. See further comments
17    in pair_mgpt.cpp and pair_mgpt.h.
18 ------------------------------------------------------------------------- */
19 
20 #ifndef READPOT__
21 #define READPOT__
22 
23 #include <cstdio>
24 #include <cstdlib>
25 #include <cstring>
26 #include "mgpt_splinetab.h"
27 
28 struct potdata {
29   double va,vb,vc,vd,ve,p1,al,rp,r00,pn;
30   double dva,dvb,dvc,dvd,dve,dp1,dal,drp,dr00;
31   double evol0,devol0;
32   double (*vpair_spline)[4],(*dvpair_spline)[4];
33   double r0,r1;
34   int nr;
35 
36   double mass,rcrit,rmax;
37 
38   int lang,lmax;
39   double anorm3,anorm4;
40   double ddl[5];
41 
42   int ipot,mode;
43   char metal[80];
44 
45   double input_vol;
46 
47   void readpot(const char *parmin_file,const char *potin_file,double vol);
48 
eval_potpotdata49   void eval_pot(double r,double *e_p,double *f_p) {
50     double d2y;
51     evalspline(nr-1,r0,r1,vpair_spline,r,e_p,f_p,&d2y);
52   }
53 
eval_virpotdata54   void eval_vir(double r,double *v_p) {
55     double dy,d2y;
56     evalspline(nr-1,r0,r1,dvpair_spline,r,v_p,&dy,&d2y);
57   }
58 };
59 
60 
61 struct potdata2 {
62   typedef double (*spline)[4];
63   spline va,vb,vc,vd,ve,p1,al,rp,r00;
64   spline dva,dvb,dvc,dvd,dve,dp1,dal,drp,dr00;
65   spline evol0,devol0;
66   double (*vpair)[4][4],(*dvpair)[4][4];
67   double r0,r1,T0,T1;
68   int nr,nt;
69 
70   potdata *potlist;
71 
72   double mass,rcrit,rmax;
73 
74   int lang,lmax;
75   spline ddl[5];
76 
77   int ipot,mode;
78   char metal[80];
79 
80   double input_vol;
81 
82   /* Functions to retrieve temperature dependent parameters */
83 
eval_tdeppotdata284   double eval_tdep(spline y,double T) {
85     double f,df,d2f;
86 
87     if(0) if(T != 3000.0)
88       printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T);
89     evalspline(nt-1,T0,T1,y,T,&f,&df,&d2f);
90     return f;
91   }
eval_tdepderivpotdata292   double eval_tdepderiv(spline y,double T) {
93     double f,df,d2f;
94 
95     if(0) if(T != 3000.0)
96       printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T);
97     evalspline(nt-1,T0,T1,y,T,&f,&df,&d2f);
98     return df;
99   }
100 
101 
102 #define make_get(param) \
103   double get_##param(double T) { return eval_tdep(param,T); } \
104   double get_d##param(double T) { return eval_tdep(d##param,T); } \
105   double get_##param##_Tderiv(double T) { return eval_tdepderiv(param,T); }
106 
107   /*
108 #define make_get(param) \
109   double get_##param(double T) { if(T != 3000.0) printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T); return potlist[3].param; } \
110   double get_d##param(double T) {  if(T != 3000.0) printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T); return potlist[3].d##param; } \
111   double get_##param##_Tderiv(double T) {  if(T != 3000.0) printf("%s:%d: Error, T = %.3f\n",__FILE__,__LINE__,T); return 0.0; }
112   */
113 
make_getpotdata2114   make_get(va)  make_get(vb)  make_get(vc)  make_get(vd)  make_get(ve)
115   make_get(p1)  make_get(al)  make_get(rp)  make_get(r00) make_get(evol0)
116 #undef make_get
117 
118   void get_anorm34(double T,double anorm3_p[1],double anorm4_p[1]) {
119     double
120       s = eval_tdep(ddl[1],T),
121       p = eval_tdep(ddl[2],T),
122       d = eval_tdep(ddl[3],T),
123       f = eval_tdep(ddl[4],T);
124     double ss = s*s, pp = p*p, dd = d*d, ff = f*f;
125     anorm3_p[0] =  s*ss + 2.0*( p*pp +  d*dd +  f*ff);
126     anorm4_p[0] = ss*ss + 2.0*(pp*pp + dd*dd + ff*ff);
127   }
get_anorm3potdata2128   double get_anorm3(double T) {
129     double a3,a4;
130     get_anorm34(T,&a3,&a4);
131     return a3;
132   }
get_anorm4potdata2133   double get_anorm4(double T) {
134     double a3,a4;
135     get_anorm34(T,&a3,&a4);
136     return a4;
137   }
get_anorm34_Tderivpotdata2138   void get_anorm34_Tderiv(double T,double danorm3_p[1],double danorm4_p[1]) {
139     double
140       s = eval_tdep(ddl[1],T),      ds = eval_tdepderiv(ddl[1],T),
141       p = eval_tdep(ddl[2],T),      dp = eval_tdepderiv(ddl[2],T),
142       d = eval_tdep(ddl[3],T),      d_d = eval_tdepderiv(ddl[3],T),
143       f = eval_tdep(ddl[4],T),      df = eval_tdepderiv(ddl[4],T);
144     double ss = s*s, pp = p*p, dd = d*d, ff = f*f;
145     danorm3_p[0] =  3.0*ds*ss + 6.0*( dp*pp +  d_d*dd +  df*ff);
146     danorm4_p[0] = 4.0*ds*s*ss + 8.0*(dp*p*pp + d_d*d*dd + df*f*ff);
147   }
get_anorm3_Tderivpotdata2148   double get_anorm3_Tderiv(double T) {
149     double da3,da4;
150     get_anorm34_Tderiv(T,&da3,&da4);
151     return da3;
152   }
get_anorm4_Tderivpotdata2153   double get_anorm4_Tderiv(double T) {
154     double da3,da4;
155     get_anorm34_Tderiv(T,&da3,&da4);
156     return da4;
157   }
158 
159   /* ... */
160 
161 
162 
parsefnamepotdata2163   char * parsefname(const char *nametemplate,int *i0,int *i1,int *stride) {
164     char *s,*p;
165 
166     if(0) {
167       s = new char[strlen(nametemplate)+1];
168     } else {
169       int len = 0;
170       while(nametemplate[len] != '\0') len = len + 1;
171       s = new char[len+1];
172     }
173     strcpy(s,nametemplate);
174 
175     p = strchr(s,'{');
176     if(p != nullptr) {
177       if(sscanf(p+1,"%d:%d:%d",i0,stride,i1) != 3) {
178         fprintf(stderr,"Error in template (\'%s\'), can not parse range.\n",nametemplate);
179         exit(1);
180       }
181       *p = '\0';
182     } else {
183       *i0 = -1;
184       *i1 = -1;
185       *stride = 1;
186     }
187     return s;
188   }
189 
maketempsplinepotdata2190   spline maketempspline(int n,potdata data[],double *ptr) {
191     int stride = &(data[1].va) - &(data[0].va);
192     spline s = new double[n-1][4];
193 
194     makespline(n,stride,ptr,s);
195     return s;
196   }
197 
readpot2potdata2198   void readpot2(const char *parmin_template,const char *potin_template,double vol) {
199     int i0,i1,stride,i0x,i1x,stridex;
200     char *parmin_file = parsefname(parmin_template,&i0 ,&i1 ,&stride );
201     char *potin_file  = parsefname( potin_template,&i0x,&i1x,&stridex);
202     int ntemp;
203 
204     potdata2 &tdeppot = *this;
205 
206     if(i0x != i0 || i1x != i1 || stridex != stride) {
207       fprintf(stderr,"Inconsistent templates. parmin_template=\'%s\', potin_template=\'%s\'\n",
208               parmin_template,potin_template);
209       exit(1);
210     }
211     if(i0 < 0 || i1 < i0 || stride <= 0 || (i1-i0)/stride+1 < 4) {
212       fprintf(stderr,"Improper temperature range. Need at least 4 temperature samples. "
213               "i0=%d,i1=%d,stride=%d,basename=\'%s\'\n",
214               i0,i1,stride,parmin_file);
215       exit(1);
216     }
217 
218     const char *parmin_suffix = strchr(parmin_template,'}')+1;
219     const char * potin_suffix = strchr( potin_template,'}')+1;
220 
221     if(parmin_suffix-1 == nullptr) {
222       fprintf(stderr,"No closing }. parmin_template=\'%s\'\n",
223               parmin_template);
224       exit(1);
225     }
226     if(potin_suffix-1 == nullptr) {
227       fprintf(stderr,"No closing }. potin_template=\'%s\'\n",
228               potin_template);
229       exit(1);
230     }
231 
232     printf("parmin_template = %s\n"
233            "parmin_file = %s\n"
234            "parmin_suffix = %s\n"
235            "T0=%d , T1=%d , stride=%d\n",
236            parmin_template,parmin_file,parmin_suffix,i0,i1,stride);
237 
238     ntemp = (i1-i0)/stride + 1;
239     /*potdata **/potlist = new potdata[ntemp];
240     char *parend = parmin_file+strlen(parmin_file);
241     char *potend = potin_file +strlen( potin_file);
242     for(int k=0; k<ntemp; k++) {
243       sprintf(parend,"%d%s",i0+k*stride,parmin_suffix);
244       sprintf(potend,"%d%s",i0+k*stride,potin_suffix);
245 
246 
247       printf("Calling readpot(%s,%s,%.3f)\n",
248              parmin_file,potin_file,vol);
249       potlist[k].readpot(parmin_file,potin_file,vol);
250 
251       if(k > 0) {
252         if(potlist[k].nr != potlist[k-1].nr) {
253           fprintf(stderr,"nr differs between file %d and %d. Exiting.\n",
254                   k,k-1);
255           exit(1);
256         }
257 
258         if(potlist[k].r0 != potlist[k-1].r0) {
259           fprintf(stderr,"r0 differs between file %d and %d. Exiting.\n",
260                   k,k-1);
261           exit(1);
262         }
263 
264         if(potlist[k].r1 != potlist[k-1].r1) {
265           fprintf(stderr,"r1 differs between file %d and %d. Exiting.\n",
266                   k,k-1);
267           exit(1);
268         }
269       }
270     }
271     tdeppot.r0 = potlist[0].r0;
272     tdeppot.r1 = potlist[0].r1;
273     tdeppot.nr = potlist[0].nr;
274     tdeppot.T0 = i0;
275     tdeppot.T1 = i1;
276     tdeppot.nt = ntemp;
277 
278     tdeppot.mass = potlist[0].mass;
279     tdeppot.rcrit = potlist[0].rcrit;
280     tdeppot.rmax = potlist[0].rmax;
281 
282     tdeppot.lang = potlist[0].lang;
283     tdeppot.lmax = potlist[0].lmax;
284     tdeppot.ipot = potlist[0].ipot;
285     tdeppot.mode = potlist[0].mode;
286     tdeppot.input_vol = potlist[0].input_vol;
287 
288     strncpy(tdeppot.metal,potlist[0].metal,sizeof(tdeppot.metal)/sizeof(char));
289     tdeppot.metal[sizeof(tdeppot.metal)/sizeof(char) - 1] = '\0';
290 
291     delete[] parmin_file;
292     delete[] potin_file;
293 
294     // Base parameters
295     tdeppot.va = maketempspline(ntemp,potlist,&(potlist[0].va));
296     tdeppot.vb = maketempspline(ntemp,potlist,&(potlist[0].vb));
297     tdeppot.vc = maketempspline(ntemp,potlist,&(potlist[0].vc));
298     tdeppot.vd = maketempspline(ntemp,potlist,&(potlist[0].vd));
299     tdeppot.ve = maketempspline(ntemp,potlist,&(potlist[0].ve));
300 
301     tdeppot.p1 = maketempspline(ntemp,potlist,&(potlist[0].p1));
302     tdeppot.al = maketempspline(ntemp,potlist,&(potlist[0].al));
303     tdeppot.rp = maketempspline(ntemp,potlist,&(potlist[0].rp));
304     tdeppot.r00 = maketempspline(ntemp,potlist,&(potlist[0].r00));
305 
306     tdeppot.evol0 = maketempspline(ntemp,potlist,&(potlist[0].evol0));
307 
308     // Volume derivatives of base parameters
309     tdeppot.dva = maketempspline(ntemp,potlist,&(potlist[0].dva));
310     tdeppot.dvb = maketempspline(ntemp,potlist,&(potlist[0].dvb));
311     tdeppot.dvc = maketempspline(ntemp,potlist,&(potlist[0].dvc));
312     tdeppot.dvd = maketempspline(ntemp,potlist,&(potlist[0].dvd));
313     tdeppot.dve = maketempspline(ntemp,potlist,&(potlist[0].dve));
314 
315     tdeppot.dp1 = maketempspline(ntemp,potlist,&(potlist[0].dp1));
316     tdeppot.dal = maketempspline(ntemp,potlist,&(potlist[0].dal));
317     tdeppot.drp = maketempspline(ntemp,potlist,&(potlist[0].drp));
318     tdeppot.dr00 = maketempspline(ntemp,potlist,&(potlist[0].dr00));
319 
320     tdeppot.devol0 = maketempspline(ntemp,potlist,&(potlist[0].devol0));
321 
322 
323     tdeppot.ddl[0] = 0;
324     for(int k = 1; k<=4; k++)
325       tdeppot.ddl[k] = maketempspline(ntemp,potlist,&(potlist[0].ddl[k]));
326 
327     {
328       double *v = new double[ntemp];
329       double (*C)[4] = new double[ntemp-1][4];
330 
331       int sz = (nr-1)*(ntemp-1);
332       //printf("Allocation:: nr=%d  ntemp=%d  size=%d\n",nr,ntemp,sz);
333       tdeppot.vpair = new double[sz][4][4];
334       tdeppot.dvpair = new double[sz][4][4];
335       /*
336       printf("vpair = %llx , dvpair = %llx",
337              (unsigned long long int) tdeppot.vpair,
338              (unsigned long long int) tdeppot.dvpair);
339       printf("   @@@@@@@@@@@@@@ nr = %d\n",nr);
340       */
341       for(int i = 0; i<nr-1; i++)
342         for(int j = 0; j<4; j++) {
343           /*
344           if(j == 5)
345             printf("    ############### i=%d\n",i);
346           */
347 
348           /* Make pair interaction interpolation functions */
349           for(int k = 0; k<ntemp; k++) {
350             if(0) if(i >= potlist[k].nr-1)
351               printf("Index error, local_nr=%d, k=%d, i=%d, nr=%d\n",nr,k,i,potlist[k].nr);
352             v[k] = potlist[k].vpair_spline[i][j];
353           }
354           makespline(ntemp,1,v,C);
355 
356           for(int k = 0; k<ntemp-1; k++)
357             for(int m = 0; m<4; m++)
358               tdeppot.vpair[k*(nr-1) + i][j][m] = C[k][m];
359 
360           /* Make pair virial interpolation functions */
361           for(int k = 0; k<ntemp; k++)
362             v[k] = potlist[k].dvpair_spline[i][j];
363           makespline(ntemp,1,v,C);
364 
365           for(int k = 0; k<ntemp-1; k++)
366             for(int m = 0; m<4; m++)
367               tdeppot.dvpair[k*(nr-1) + i][j][m] = C[k][m];
368 
369         }
370 
371       delete[] C;
372       delete[] v;
373     }
374 
375   }
376 
eval2Dsplinepotdata2377   void eval2Dspline(double fun[][4][4],double r,double T,double *e_p,double *f_p,double *dedT_p) {
378     double C[4],dC[4],dd;
379     double Tfrac = (T-T0)/(T1-T0) * (nt-1);
380     double rfrac = (r-r0)/(r1-r0) * (nr-1);
381     int k = (int) Tfrac;
382     int i = (int) rfrac;
383     int j;
384 
385     if(k < 0) k = 0;
386     else if(k > nt-2) k = nt-2;
387 
388     if(i < 0) i = 0;
389     else if(i > nr-2) i = nr-2;
390 
391     /*
392     printf("eval_pot nr=%d  nt=%d\n",nr,nt);
393     printf("eval_pot r=%.3f  T=%.3f  k=%d  i=%d\n",
394            r,T,k,i);
395     printf("Tfrac=%.3f  Tfrac-k=%.3f\n",Tfrac,Tfrac-k);
396     printf("rfrac=%.3f  rfrac-i=%.3f\n",rfrac,rfrac-i);
397     */
398     for(j = 0; j<4; j++) {
399       evalcubic(fun[k*(nr-1) + i][j],Tfrac-k,&C[j],&dC[j],&dd);
400       dC[j] = dC[j] * ((nt-1) / (T1-T0));
401     }
402     /*
403     printf("C coeff: %.3e  %.3e  %.3e  %.3e\n",
404            C[0],C[1],C[2],C[3]);
405     */
406     evalcubic(C,rfrac-i,e_p,f_p,&dd);
407     evalcubic(dC,rfrac-i,dedT_p,&dd,&dd);
408     *f_p *= (nr-1) / (r1-r0);
409   }
eval_potpotdata2410   void eval_pot(double r,double T,double *e_p,double *f_p,double *dedT_p) {
411     eval2Dspline(vpair,r,T,e_p,f_p,dedT_p);
412   }
eval_virpotdata2413   void eval_vir(double r,double T,double *v_p) {
414     double vf,dvdT;
415     eval2Dspline(dvpair,r,T,v_p,&vf,&dvdT);
416   }
417 };
418 
419 
420 #endif
421