Commit JT 042519

- moved examples/SPIN/gneb/interpolate/ toward
tools/spin/interpolate/gneb
- added gneb in the command files in doc/src
- modified the error messages in neb_spin.h/cpp
This commit is contained in:
julient31
2019-04-25 09:18:33 -06:00
parent b8ae466999
commit 7c3f619e13
13 changed files with 19 additions and 12 deletions

View File

@ -1,31 +0,0 @@
Interpolate.x tries to perform a cubic polynomial interpolation
of the MEP.
Compile the program with:
gcc interpolate.c -o interpolate.x -lm -lgsl
and then run it as:
./interpolate.x
Enter N-1 (with N the number of replica).
The program reads the "neb_init.dat", in which you need to
replace the current values by your GNEB outputs.
Each line corresponds to a replica, and has to respectively
contain:
Reac. coords /tab/ pot. energy /tab/ fm dot tangent /tab/
geodesic dist to next replica
All those information can be provided by the verbose output of
a neb/spin calculation
The progam outputs the interpolation result, and the
interpolated MEP in "interpolation_result.dat".
This code is a courtesy of Aleksei Ivanov, University of
Iceland.
For more explanation about this calculation, see Appendix D
of the following reference:
Bessarab, P. F., Uzdin, V. M., & Jónsson, H. (2015).
Computer Physics Communications, 196, 335-347.

View File

@ -1,108 +0,0 @@
/* ------------------------------------------------------------------------
Provide some explanation here
------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------
This program is a courtesy of Aleksei Ivanov (Univ. of Iceland)
Contributing authors: Aleksei Ivanov (Univ. of Iceland),
Julien Tranchida (SNL)
------------------------------------------------------------------------- */
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <stdlib.h>
// calculate cubic coefficients
void count_coefficient(double *V, double *F, double *R, double *a, double *b,
double *c,double *d,int M){
/* R = square of distance between images*/
/* V = energy of images */
/* F = projection of real real forces along the path? */
int i;
for(i = 0; i < M ; i++) {
a[i] = ( -2.0*(V[i+1]-V[i])/R[i] - F[i] - F[i+1] ) / (R[i]*R[i]);
//a[i] = ( -2.0*(V[i+1]-V[i])/sqrt(R[i]) - F[i] - F[i+1] ) / R[i];
b[i] = ( 3.0*(V[i+1]-V[i]) + (2.0*F[i]+F[i+1])*R[i] ) / (R[i]*R[i]);
//b[i] = ( 3.0*(V[i+1]-V[i]) + (2.0*F[i]+F[i+1])*sqrt(R[i]) ) / R[i];
c[i] = -F[i];
d[i] = V[i];
}
}
// cubic spline
double spline(double a,double b,double c,double d,double x) {
return (a*x*x*x + b*x*x + c*x + d);
}
int main() {
int M=0; // M+1 = number of images
double *fmdottan; // projection of real forces on tangent path
double *coords; // initial value of reaction coords
double *V; // energy of images
double *dist; // square of the distance between images
double *a, *b, *c, *d ; // coefficients of cubic functions
double x; // reaction coordinate
double A,B; // additional variables for rnd
double length = 0.0;
int i,j;
FILE *data;
printf("Enter M = number of images - 1 \n");
scanf("%d",&M);
// allocating and initializing tables
a = calloc(M,sizeof(double)); // cubic coefficients
b = calloc(M,sizeof(double));
c = calloc(M,sizeof(double));
d = calloc(M,sizeof(double));
V = calloc((M+1),sizeof(double)); // energies
coords = calloc((M+1),sizeof(double)); // reaction coordinates
fmdottan = calloc((M+1),sizeof(double)); // fm dot tangent
dist = calloc(M+1,sizeof(double)); // distance between images
// reading input file
if((data=fopen("neb_init.dat","r")) == NULL) {
printf("Incorrect input file name.");
return 0;
}
for(j=0; j < M+1; j++) {
fscanf(data,"%lf\t%lf\t%lf\t%lf\n",&coords[j],&V[j],&fmdottan[j],&dist[j]);
length += dist[j];
printf("%lf %lf %lf %lf\n",coords[j],V[j],fmdottan[j],dist[j]);
}
if( (fclose(data)) == 0) {
printf("Data stored, input file closed.\n ");
}
// calculate value of coefficients
count_coefficient(V,fmdottan,dist,a,b,c,d,M);
// plot result of the interpolation
if( ( data=fopen("interpolation_result.dat","w") )== NULL) {
printf("Interpolation file could not be open.");
return 0;
}
A = B = 0.0;
for(i = 0; i < M ; i++) {
B += dist[i];
printf("%13le\n",B);
for(j = 0; j <= 1000; j++) {
x = dist[i]*1.0e-3*j;
fprintf(data,"%13lf\t%13le\n",(x+A)/length,spline(a[i],b[i],c[i],d[i],x));
}
A += dist[i];
}
return 0;
}

View File

@ -1,8 +0,0 @@
0.0000000 -3.915271 3.4995081e-17 2.4573077
0.14285714 -3.9148148 -0.00059075739 2.4573077
0.28571429 -3.9136926 -0.00072315767 2.4573077
0.42857143 -3.9127883 -0.0003191228 2.4573265
0.57142857 -3.9127883 0.0003191228 2.4403341
0.71428571 -3.9136926 0.00072315767 2.4044093
0.85714286 -3.9148148 0.00059075739 2.3766041
1.0000000 -3.915271 -4.1231828e-17 0.0000000