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