Merge pull request #2258 from akohlmey/phana-update (2923dcbb) · Commits · 郑智淋 / lammps

tools/phonon/Makefile

+26 −32

Original line number	Original line	Diff line number	Diff line
	.SUFFIXES : .o .cpp		.SUFFIXES : .o .cpp
	# compiler and flags		# compiler and flags
	CC = g++ -Wall		CC = g++ -Wno-unused-result
	LINK = $(CC)		LINK = $(CC) -static
	CFLAGS = -O3 $(DEBUG) $(UFLAG)		CFLAGS = -O3 $(DEBUG) $(UFLAG)
	#
	OFLAGS = -O3 $(DEBUG)		OFLAGS = -O3 $(DEBUG)
	INC = $(LPKINC) $(TCINC) $(SPGINC)		INC = $(LPKINC) $(TCINC) $(SPGINC) $(FFTINC)
	LIB = $(LPKLIB) $(TCLIB) $(SPGLIB)		LIB = $(LPKLIB) $(TCLIB) $(SPGLIB) $(FFTLIB)
	#
	# cLapack library needed		# cLapack library needed
	LPKINC =		LPKINC = -I/opt/clapack/3.2.1/include
	LPKLIB =-llapack		LPKLIB = -L/opt/clapack/3.2.1/lib -lclapack -lblas -lf2c -lm
	#
	#		# Tricubic library needed
	# spglib 1.8.2, used to get the irreducible q-points		TCINC = -I/opt/tricubic/1.0/include
	# if UFLAG is not set, spglib won't be used.		TCLIB = -L/opt/tricubic/1.0/lib -ltricubic

	# UFLAG = -DUseSPG		# spglib, used to get the irreducible q-points
	# SPGINC = -I/opt/libs/spglib/1.8.2/include		# if SFLAG is not set, spglib won't be used.
	# SPGLIB = -L/opt/libs/spglib/1.8.2/lib -lsymspg		SFLAG = -DUseSPG
			SPGINC = -I/opt/spglib/1.9.7/include/spglib
			SPGLIB = -L/opt/spglib/1.9.7/lib -lsymspg

	# if spglib other than version 1.8.2 is used, please		# FFTW 3， used to deduce the force constants in real space
	# modify file phonon.cpp, instruction can be found by searching 1.8.2		# if FFLAG is not set, fftw won't be used.
			FFLAG = -DFFTW3
			FFTINC = -I/opt/fftw/3.3.7/include
			FFTLIB = -L/opt/fftw/3.3.7/lib -lfftw3

	# Debug flags		# Debug flags
	# DEBUG = -g -DDEBUG		# DEBUG = -g -DDEBUG
			UFLAG = $(SFLAG) $(FFLAG)

	#====================================================================		#====================================================================
	ROOT = phana		ROOT = phana
	# executable name		# executable name
	@@ -44,7 +51,7 @@ clean:
	rm -f .o ~ *.mod ${EXE}		rm -f .o ~ *.mod ${EXE}

	tar:		tar:
	rm -f ${ROOT}.tar; tar -czvf ${ROOT}.tar.gz .cpp .h Makefile README		rm -f ${ROOT}.tar.gz; tar -czvf ${ROOT}.tar.gz .cpp .h Makefile README

	ver:		ver:
	@echo "#define VERSION `git log\|grep '^commit'\|wc -l`" > version.h		@echo "#define VERSION `git log\|grep '^commit'\|wc -l`" > version.h
	@@ -58,16 +65,3 @@ ver:
	$(CC) $(CFLAGS) -c $<		$(CC) $(CFLAGS) -c $<
	.cpp.o:		.cpp.o:
	$(CC) $(CFLAGS) $(INC) -c $<		$(CC) $(CFLAGS) $(INC) -c $<

	#====================================================================
	# dependencies
	disp.o: disp.cpp phonon.h dynmat.h memory.h interpolate.h green.h timer.h \
	global.h
	dynmat.o: dynmat.cpp dynmat.h memory.h interpolate.h version.h global.h
	green.o: green.cpp green.h memory.h global.h
	interpolate.o: interpolate.cpp interpolate.h memory.h global.h
	main.o: main.cpp dynmat.h memory.h interpolate.h phonon.h
	memory.o: memory.cpp memory.h
	phonon.o: phonon.cpp phonon.h dynmat.h memory.h interpolate.h green.h \
	timer.h global.h
	timer.o: timer.cpp timer.h

tools/phonon/README

+18 −9

Original line number	Original line	Diff line number	Diff line
	@@ -5,9 +5,13 @@
	analyse the phonon related information.		analyse the phonon related information.
	#-------------------------------------------------------------------------------		#-------------------------------------------------------------------------------
	1. Dependencies		1. Dependencies
	The LAPACK library is needed to solve the eigen problems.		The clapack library is needed to solve the eigen problems,
	http://www.netlib.org/lapack/		which could be downloaded from:
	Intel MKL can be used as well.		http://www.netlib.org/clapack/

			The tricubic library is also needed to do tricubic interpolations,
			which could now be obtained from:
			https://github.com/nbigaouette/libtricubic/

	The spglib is optionally needed, enabling one to evaluate the		The spglib is optionally needed, enabling one to evaluate the
	phonon density of states or vibrational thermal properties		phonon density of states or vibrational thermal properties
	@@ -17,6 +21,12 @@
	It can be obtained from:		It can be obtained from:
	http://spglib.sourceforge.net/		http://spglib.sourceforge.net/

			FFTW 3 might also be needed if you would like to interface with
			phonopy: necessary input files for phonopy will be prepared so
			that you can make use of the functions provided by phonopy.
			FFTW 3 can be downloaded from:
			http://www.fftw.org

	2. Compilation		2. Compilation
	To compile the code, one needs therefore to install the above		To compile the code, one needs therefore to install the above
	libraries and set the paths correctly in the Makefile.		libraries and set the paths correctly in the Makefile.
	@@ -26,17 +36,16 @@

	3. Unit system		3. Unit system
	The units of the output frequencies by this code is THz for		The units of the output frequencies by this code is THz for
	LAMMPS units "real", "si", "metal", and "cgs"; in these cases,		LAMMPS units "real", "si", "metal", "cgs", "micro", "nano";
	the frequencies are $\nu$ instead of $\omega$.		in these cases, the frequencies are $\nu$ instead of $\omega$.

	4. Updates		4. Updates
	For updates of phana, please check:		For updates of phana, please check:
	http://nes.sjtu.edu.cn/english/research/software.htm
	or
	https://github.com/lingtikong/phana.git		https://github.com/lingtikong/phana.git

	5. Bug report		5. Bug report
	If any bug found, please drop a line to: konglt(at)sjtu.edu.cn		If any bug found, please drop a line to: konglt(at)sjtu.edu.cn

	#-------------------------------------------------------------------------------		#-------------------------------------------------------------------------------
	Author: Ling-Ti Kong, konglt(at)sjtu.edu.cn		Author: Ling-Ti Kong, konglt(at)sjtu.edu.cn
	Oct 2015		May 2020

tools/phonon/disp.cpp

+135 −2832

File changed.

Preview size limit exceeded, changes collapsed.

tools/phonon/dynmat.cpp

+568 −473

Original line number	Original line	Diff line number	Diff line
	@@ -3,12 +3,10 @@
	#include "version.h"		#include "version.h"
	#include "global.h"		#include "global.h"

	extern "C" void zheevd_(char , char , long int , doublecomplex ,		/* ----------------------------------------------------------------------------
	long int , double , doublecomplex *,		* Class DynMat stores the Dynamic Matrix read from the binary file from
	long int , double , long int , long int ,		* fix-phonon and incorporates some operations w.r.t the DM.
	long int , long int );		* ---------------------------------------------------------------------------- */

	// to initialize the class
	DynMat::DynMat(int narg, char **arg)		DynMat::DynMat(int narg, char **arg)
	{		{
	attyp = NULL;		attyp = NULL;
	@@ -69,42 +67,48 @@ DynMat::DynMat(int narg, char **arg)
	}		}

	// read header info from the binary file		// read header info from the binary file
	if ( fread(&sysdim, sizeof(int), 1, fp) != 1) {printf("\nError while reading sysdim from file: %s\n", binfile); fclose(fp); exit(2);}		if (fread(&sysdim, sizeof(int), 1, fp) != 1){
	if ( fread(&nx, sizeof(int), 1, fp) != 1) {printf("\nError while reading nx from file: %s\n", binfile); fclose(fp); exit(2);}		printf("\nError while reading sysdim from file: %s\n", binfile);
	if ( fread(&ny, sizeof(int), 1, fp) != 1) {printf("\nError while reading ny from file: %s\n", binfile); fclose(fp); exit(2);}		fclose(fp);
	if ( fread(&nz, sizeof(int), 1, fp) != 1) {printf("\nError while reading nz from file: %s\n", binfile); fclose(fp); exit(2);}		exit(2);
	if ( fread(&nucell, sizeof(int), 1, fp) != 1) {printf("\nError while reading nucell from file: %s\n", binfile); fclose(fp); exit(2);}		}
	if ( fread(&boltz, sizeof(double), 1, fp) != 1) {printf("\nError while reading boltz from file: %s\n", binfile); fclose(fp); exit(2);}		if (fread(&nx, sizeof(int), 1, fp) != 1){
			printf("\nError while reading nx from file: %s\n", binfile);
	fftdim = sysdimnucell; fftdim2 = fftdimfftdim;		fclose(fp);
			exit(2);
			}
			if (fread(&ny, sizeof(int), 1, fp) != 1){
			printf("\nError while reading ny from file: %s\n", binfile);
			fclose(fp);
			exit(2);
			}
			if (fread(&nz, sizeof(int), 1, fp) != 1){
			printf("\nError while reading nz from file: %s\n", binfile);
			fclose(fp);
			exit(2);
			}
			if (fread(&nucell, sizeof(int), 1, fp) != 1){
			printf("\nError while reading nucell from file: %s\n", binfile);
			fclose(fp);
			exit(2);
			}
			if (fread(&boltz, sizeof(double), 1, fp) != 1){
			printf("\nError while reading boltz from file: %s\n", binfile);
			fclose(fp);
			exit(2);
			}

			fftdim = sysdim * nucell;
			fftdim2 = fftdim * fftdim;
	npt = nx * ny * nz;		npt = nx * ny * nz;

	// display info related to the read file		// display info related to the read file
	printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");		ShowInfo();
	printf("Dynamical matrix is read from file: %s\n", binfile);
	printf("The system size in three dimension: %d x %d x %d\n", nx, ny, nz);
	printf("Number of atoms per unit cell : %d\n", nucell);
	printf("System dimension : %d\n", sysdim);
	printf("Boltzmann constant in used units : %g\n", boltz);
	for (int i = 0; i < 80; ++i) printf("=");
	printf("\n");
	if (sysdim < 1\|\|sysdim > 3\|\|nx < 1\|\|ny < 1\|\|nz < 1\|\|nucell < 1){		if (sysdim < 1\|\|sysdim > 3\|\|nx < 1\|\|ny < 1\|\|nz < 1\|\|nucell < 1){
	printf("Wrong values read from header of file: %s, please check the binary file!\n", binfile);		printf("Wrong values read from header of file: %s, please check the binary file!\n", binfile);
	fclose(fp); exit(3);		fclose(fp); exit(3);
	}		}
			Define_Conversion_Factor();
	funit = new char[4];
	strcpy(funit, "THz");
	if (boltz == 1.){eml2f = 1.; delete funit; funit = new char[27]; strcpy(funit,"sqrt(epsilon/(m.sigma^2))");}
	else if (boltz == 0.0019872067) eml2f = 3.256576161;
	else if (boltz == 8.617343e-5) eml2f = 15.63312493;
	else if (boltz == 1.3806504e-23) eml2f = 1.;
	else if (boltz == 1.3806504e-16) eml2f = 1.591549431e-14;
	else {
	printf("WARNING: Because of float precision, I cannot get the factor to convert sqrt(E/ML^2)\n");
	printf("into THz, instead, I set it to be 1; you should check the unit used by LAMMPS.\n");
	eml2f = 1.;
	}

	// now to allocate memory for DM		// now to allocate memory for DM
	memory = new Memory();		memory = new Memory();
	@@ -123,11 +127,31 @@ DynMat::DynMat(int narg, char **arg)
	memory->create(attyp, nucell, "DynMat:attyp");		memory->create(attyp, nucell, "DynMat:attyp");
	memory->create(M_inv_sqrt, nucell, "DynMat:M_inv_sqrt");		memory->create(M_inv_sqrt, nucell, "DynMat:M_inv_sqrt");

	if ( (int) fread(&Tmeasure, sizeof(double), 1, fp) != 1 ){printf("\nError while reading temperature from file: %s\n", binfile); fclose(fp); exit(3);}		if (fread(&Tmeasure, sizeof(double), 1, fp) != 1){
	if ( (int) fread(&basevec[0], sizeof(double), 9, fp) != 9 ){printf("\nError while reading lattice info from file: %s\n", binfile); fclose(fp); exit(3);}		printf("\nError while reading temperature from file: %s\n", binfile);
	if ( (int) fread(basis[0], sizeof(double), fftdim, fp) != fftdim){printf("\nError while reading basis info from file: %s\n", binfile); fclose(fp); exit(3);}		fclose(fp);
	if ( (int) fread(&attyp[0], sizeof(int), nucell, fp) != nucell){printf("\nError while reading atom types from file: %s\n", binfile); fclose(fp); exit(3);}		exit(3);
	if ( (int) fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != nucell){printf("\nError while reading atomic masses from file: %s\n", binfile); fclose(fp); exit(3);}		}
			if (fread(&basevec[0], sizeof(double), 9, fp) != 9){
			printf("\nError while reading lattice info from file: %s\n", binfile);
			fclose(fp);
			exit(3);
			}
			if (fread(basis[0], sizeof(double), fftdim, fp) != fftdim){
			printf("\nError while reading basis info from file: %s\n", binfile);
			fclose(fp);
			exit(3);
			}
			if (fread(&attyp[0], sizeof(int), nucell, fp) != nucell){
			printf("\nError while reading atom types from file: %s\n", binfile);
			fclose(fp);
			exit(3);
			}
			if (fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != nucell){
			printf("\nError while reading atomic masses from file: %s\n", binfile);
			fclose(fp);
			exit(3);
			}
	fclose(fp);		fclose(fp);

	car2dir();		car2dir();
	@@ -158,7 +182,10 @@ DynMat::DynMat(int narg, char **arg)
	return;		return;
	}		}

	// to destroy the class
			/* ----------------------------------------------------------------------------
			* Free the memories used.
			* ---------------------------------------------------------------------------- */
	DynMat::~DynMat()		DynMat::~DynMat()
	{		{
	// destroy all memory allocated		// destroy all memory allocated
	@@ -172,7 +199,9 @@ DynMat::~DynMat()
	memory->destroy(basis);		memory->destroy(basis);
	memory->destroy(DM_all);		memory->destroy(DM_all);
	memory->destroy(M_inv_sqrt);		memory->destroy(M_inv_sqrt);
	if (memory) delete memory;		delete memory;

			return;
	}		}

	/* ----------------------------------------------------------------------------		/* ----------------------------------------------------------------------------
	@@ -209,6 +238,7 @@ void DynMat::writeDMq(double *q)
	}		}
	fprintf(fp,"\n");		fprintf(fp,"\n");
	fclose(fp);		fclose(fp);

	return;		return;
	}		}

	@@ -217,7 +247,6 @@ return;
	* ---------------------------------------------------------------------------- */		* ---------------------------------------------------------------------------- */
	void DynMat::writeDMq(double q, const double qr, FILE fp)		void DynMat::writeDMq(double q, const double qr, FILE fp)
	{		{

	fprintf(fp, "%lg %lg %lg %lg ", q[0], q[1], q[2], qr);		fprintf(fp, "%lg %lg %lg %lg ", q[0], q[1], q[2], qr);

	for (int i = 0; i < fftdim; ++i)		for (int i = 0; i < fftdim; ++i)
	@@ -235,9 +264,9 @@ return;
	int DynMat::geteigen(double *egv, int flag)		int DynMat::geteigen(double *egv, int flag)
	{		{
	char jobz, uplo;		char jobz, uplo;
	long int n, lda, lwork, lrwork, *iwork, liwork, info;		integer n, lda, lwork, lrwork, *iwork, liwork, info;
	doublecomplex *work;		doublecomplex *work;
	double w = &egv[0], rwork;		doublereal w = &egv[0], rwork;

	n = fftdim;		n = fftdim;
	if (flag) jobz = 'V';		if (flag) jobz = 'V';
	@@ -276,6 +305,7 @@ return info;
	void DynMat::getDMq(double *q)		void DynMat::getDMq(double *q)
	{		{
	interpolate->execute(q, DM_q[0]);		interpolate->execute(q, DM_q[0]);

	return;		return;
	}		}

	@@ -344,8 +374,7 @@ void DynMat::EnforceASR()
	char *ptr = strtok(str," \t\n\r\f");		char *ptr = strtok(str," \t\n\r\f");
	if (ptr) nasr = atoi(ptr);		if (ptr) nasr = atoi(ptr);
	if (nasr < 1){		if (nasr < 1){
	for (int i=0; i<80; i++) printf("=");		for (int i=0; i<80; i++) printf("="); printf("\n");
	printf("\n");
	return;		return;
	}		}

	@@ -408,11 +437,10 @@ void DynMat::EnforceASR()
	for (int i = 0; i < fftdim; ++i){		for (int i = 0; i < fftdim; ++i){
	printf("%lg ", egvs[i]);		printf("%lg ", egvs[i]);
	if (i%10 == 9) printf("\n");		if (i%10 == 9) printf("\n");
	if (i == 99){ printf("...... (%d more skipped)", fftdim-100); break;}		if (i == 99){ printf("...... (%d more skiped)", fftdim-100); break;}
	}		}
	printf("\n");		printf("\n");
	for (int i = 0; i < 80; ++i) printf("=");		for (int i = 0; i < 80; ++i) printf("="); printf("\n\n");
	printf("\n\n");

	return;		return;
	}		}
	@@ -464,7 +492,7 @@ return;
	* --------------------------------------------------------------------*/		* --------------------------------------------------------------------*/
	void DynMat::GaussJordan(int n, double *Mat)		void DynMat::GaussJordan(int n, double *Mat)
	{		{
	int i,icol=0,irow=0,j,k,l,ll,idr,idc;		int i,icol,irow,j,k,l,ll,idr,idc;
	int indxc,indxr,*ipiv;		int indxc,indxr,*ipiv;
	double big, nmjk;		double big, nmjk;
	double dum, pivinv;		double dum, pivinv;
	@@ -591,7 +619,74 @@ void DynMat::ShowVersion()
	printf(" )___/ ) _ ( /(__)\\ ) ( /(__)\\ \n");		printf(" )___/ ) _ ( /(__)\\ ) ( /(__)\\ \n");
	printf(" (__) (_) (_)(__)(__)(_)\\_)(__)(__)\n");		printf(" (__) (_) (_)(__)(__)(_)\\_)(__)(__)\n");
	printf("\nPHonon ANAlyzer for Fix-Phonon, version 2.%02d, compiled on %s.\n", VERSION, __DATE__);		printf("\nPHonon ANAlyzer for Fix-Phonon, version 2.%02d, compiled on %s.\n", VERSION, __DATE__);
			printf("Reference: https://doi.org/10.1016/j.cpc.2011.04.019\n");

			return;
			}

			/* ----------------------------------------------------------------------------
			* Private method to define the conversion factor from the DM measured to THZ
			* for the phonon frequency and force constants.
			* ---------------------------------------------------------------------------- */
			void DynMat::Define_Conversion_Factor()
			{
			funit = new char[4];
			strcpy(funit, "THz");

			if (fabs(boltz - 1.) <= ZERO){ // LJ Unit
			eml2f = eml2fc = 1.;
			delete funit;
			funit = new char[27];
			strcpy(funit,"sqrt(epsilon/(m.sigma^2))");

			} else if (fabs(boltz - 0.0019872067) <= ZERO){ // real
			eml2f = 3.255487031;
			eml2fc = 0.0433641042418;

			} else if (fabs(boltz*1.e3 - 8.617343e-2) <= ZERO){ // metal
			eml2f = 15.633304237154924;
			eml2fc = 1.;

			} else if (fabs(boltz*1.e20 - 1.3806504e-3) <= ZERO){ // si
			eml2f = 1.591549431e-13;
			eml2fc = 0.06241509074460763;

			} else if (fabs(boltz*1.e13 - 1.3806504e-3) <= ZERO){ // cgs
			eml2f = 1.591549431e-13;
			eml2fc = 6.241509074460763e-05;

			} else if (fabs(boltz*1.e3 - 3.16681534e-3) <= ZERO){ // electron
			eml2f = 154.10792761319672;
			eml2fc = 97.1736242922823;

			} else if (fabs(boltz*1.e5 - 1.3806504e-3) <= ZERO){ // micro
			eml2f = 1.5915494309189532e-07;
			eml2fc = 6.241509074460763e-05;

			} else if (fabs(boltz - 0.013806504) <= ZERO){ // nano
			eml2f = 0.0001591549431;
			eml2fc = 6.241509074460763e-05;

			} else {
			printf("WARNING: Perhaps because of float precision, I cannot get the factor to convert\n");
			printf("sqrt(E/ML^2)/(2*pi) into THz, instead, I set it to 1; you should check the unit\nused by LAMMPS.\n");
			eml2f = eml2fc = 1.;
			}
			return;
			}

			/* ----------------------------------------------------------------------------
			* Private method to output the information read
			* ---------------------------------------------------------------------------- */
			void DynMat::ShowInfo()
			{
			printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");
			printf("Dynamical matrix is read from file: %s\n", binfile);
			printf("The system size in three dimension: %d x %d x %d\n", nx, ny, nz);
			printf("Number of atoms per unit cell : %d\n", nucell);
			printf("System dimension : %d\n", sysdim);
			printf("Boltzmann constant in used units : %g\n", boltz);
			for (int i = 0; i < 80; ++i) printf("="); printf("\n");
	return;		return;
	}		}
	/* --------------------------------------------------------------------*/		/* --------------------------------------------------------------------*/

tools/phonon/dynmat.h

+5 −5

Original line number	Original line	Diff line number	Diff line
	@@ -7,8 +7,6 @@
	#include "memory.h"		#include "memory.h"
	#include "interpolate.h"		#include "interpolate.h"

	using namespace std;

	class DynMat {		class DynMat {
	public:		public:

	@@ -17,7 +15,7 @@ public:

	int nx, ny, nz, nucell;		int nx, ny, nz, nucell;
	int sysdim, fftdim;		int sysdim, fftdim;
	double eml2f;		double eml2f, eml2fc;
	char *funit;		char *funit;

	void getDMq(double *);		void getDMq(double *);
	@@ -30,7 +28,9 @@ public:
	doublecomplex **DM_q;		doublecomplex **DM_q;

	int flag_latinfo;		int flag_latinfo;
			int npt, fftdim2;
	double Tmeasure, basevec[9], ibasevec[9];		double Tmeasure, basevec[9], ibasevec[9];
			double *M_inv_sqrt;
	double **basis;		double **basis;
	int *attyp;		int *attyp;

	@@ -40,14 +40,12 @@ private:
	Interpolate *interpolate;		Interpolate *interpolate;

	Memory *memory;		Memory *memory;
	int npt, fftdim2;

	int nasr;		int nasr;
	void EnforceASR();		void EnforceASR();

	char binfile, dmfile;		char binfile, dmfile;
	double boltz, q[3];		double boltz, q[3];
	double *M_inv_sqrt;

	doublecomplex **DM_all;		doublecomplex **DM_all;

	@@ -56,6 +54,8 @@ private:
	void GaussJordan(int, double *);		void GaussJordan(int, double *);

	void help();		void help();
			void ShowInfo();
	void ShowVersion();		void ShowVersion();
			void Define_Conversion_Factor();
	};		};
	#endif		#endif

Admin message