/* Interpolates 2 POSCAR files for elastic band method. */

/* Edit the number below to the number of different atom types present. */

#define ATOMTYPES 4
   
/*
   Program assumes 2 files are concatinated, with an aditional
   line at the beginning of the file. This line contains the number
   of subdirectories to be written to, followed by a 'T' or
   an 'F', depending on whether or not the center of mass is shifted,
   followed by another 'T' or 'F' depending on whether or not the
   unit cell axes are interpolated. 
   If the word 'Selective' is chosen, then 'dynamics' must also follow.
   The program creates the required subdirectories automatically, and
   also warns if coordinate errors are perhaps present due to the wrap 
   around periodicity in the unit cell.
   To run, type the program name followed by the concatinated input file.  
   Scott R. Shannon, 14/10/98 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#define CUTOFF 0.5
#define MAXATOM 200

main(int argc, char *argv[])
{
	int i,j,k,dum,total,select=0,reps,numatoms[ATOMTYPES];
	FILE *indat,*outdat,*ifp;
	double latc,uaxis[2][3][3],coords[2][MAXATOM][3];
	double increments[MAXATOM][3]={0.0},axisinc[3][3]={0.0};
	double cm[2][3]={0.0},cmshift[3]={0.0};
	char string[50],doshift[50],interuaxis[50],command[50];
	char title[50],cotype[50],filename[50];
	char selstring[MAXATOM][3][5];
	
	if ( (indat = fopen(argv[1],"r"))==NULL) { printf("\nIncorrect data file name\n\n");
						exit(); }
	fscanf(indat,"%d %s %s",&reps,doshift,interuaxis);

	for (i=0;i<2;++i) {

	fscanf(indat,"%s",title);  /* title */

	fscanf(indat,"%lf",&latc);  /* lattice constants */

	for (j=0;j<3;++j)
		for (k=0;k<3;++k)
			fscanf(indat,"%lf",&uaxis[i][j][k]); /* unit vectors */

	total = 0;
	for (j=0;j<ATOMTYPES;++j) {
		fscanf(indat,"%d",&numatoms[j]);  /* number of each atoms */
		total += numatoms[j]; }

	fscanf(indat,"%s",string);  /* read in type */

	select = FALSE;
	if (strcmp(string,"Selective")==0) {
			select = TRUE;
			fscanf(indat,"%s",string); /* read in "Dynamics" */
			fscanf(indat,"%s",cotype);  /* read in "Direct|Cartesian" */
			}

	for (j=0;j<total;++j) {
		for (k=0;k<3;++k) {
			if (fscanf(indat,"%lf",&coords[i][j][k])==EOF) {
				printf("\n\nError: end of file reached\n");
				printf("Check number of atoms. Exiting\n\n");
				exit(); }
			cm[i][k] += coords[i][j][k];
			}
		if (select)
			fscanf(indat,"%s %s %s",selstring[j][0],
					     selstring[j][1],
					     selstring[j][2]);
		}
		
	}
	
/*	for (i=0;i<total;++i)
		printf("%f %f %f\n",coords[1][i][0],coords[1][i][1],
						coords[1][i][2]);
*/

	for (i=0;i<total;++i) /* check for wrap around errors in POSCAR's */
	   for (j=0;j<3;++j)
		if (fabs(coords[0][i][j] - coords[1][i][j]) > CUTOFF){
			printf("\n\nWarning! Possible error due to wrap around in unit cell\n");
			printf("\nCheck atom %d\n\n",i+1);
			}
		
	if (strcmp(doshift,"T")==0) {  /* shift center of mass */
		for (i=0;i<3;++i)
			cmshift[i] = cm[0][i]/total - cm[1][i]/total;
		printf("\nCenter of mass for second cell will be shifted by:  %e  %e  %e\n\n",
			cmshift[0],cmshift[1],cmshift[2]);
		for (i=0;i<3;++i)
			cmshift[i] = cmshift[i]/(reps-1);
		}

	if (strcmp(interuaxis,"T")==0) {  /* increment uc axes */
		for (i=0;i<3;++i)
		for (j=0;j<3;++j) 
			axisinc[i][j] =
			  (uaxis[1][i][j] - uaxis[0][i][j])/(reps-1);
		}

	for (i=0;i<total;++i)         /* calculate coord increments */
		for (j=0;j<3;++j) 
		     increments[i][j] = 
			(coords[1][i][j] - coords[0][i][j])/(reps-1);

	for (i=0;i<reps;++i) {   /* for each directory.. */
	
	sprintf(command,"find . -name 0%d",i);
	ifp = popen(command,"r");
	if (getc(ifp)==EOF) {     /* if directory not present, create it */
		sprintf(command,"mkdir 0%d",i);
		system(command);  } 

	sprintf(filename,"./0%d/POSCAR",i);
	outdat = fopen(filename,"w");
	
	fprintf(outdat,"%s\n%f\n",title,latc);
	
	for (j=0;j<3;++j) {   /* output new unit cell axes */
	       	 for (k=0;k<3;++k)
		 	fprintf(outdat,"%f  ",uaxis[0][j][k] + i*axisinc[j][k]);
		 fprintf(outdat,"\n");
		 }

	for (j=0;j<ATOMTYPES;++j)
		fprintf(outdat,"%d  ",numatoms[j]);
	fprintf(outdat,"\n");
	
	if (select) fprintf(outdat,"Selective Dynamics\n");
	
	fprintf(outdat,"%s\n",cotype);

	for (j=0;j<total;++j) {  /* print out coords and selective ops */
	
	for (k=0;k<3;++k) 
		fprintf(outdat,"%f  ",
			coords[0][j][k] + i*cmshift[k] + i*increments[j][k]);

	if (select) {
		for (k=0;k<3;++k)
			fprintf(outdat,"%s ",selstring[j][k]); }
	fprintf(outdat,"\n");	

	}  /* end coords */
	
	fclose(outdat);
	pclose(ifp);
	
 	}		
					
}