Update DihedralTable (91e84693) · Commits · 郑智淋 / lammps

src/USER-MISC/dihedral_table.cpp

+61 −99

Original line number	Diff line number	Diff line
		@@ -36,6 +36,9 @@
		#include "utils.h"
		#include "dihedral_table.h"
		#include "utils.h"
		#include "tokenizer.h"
		#include "table_file_reader.h"
		#include "fmt/format.h"

		#include "math_const.h"
		#include "math_extra.h"
		@@ -458,8 +461,7 @@ static double Phi(double const *x1, //array holding x,y,z coords atom 1
		DihedralTable::DihedralTable(LAMMPS *lmp) : Dihedral(lmp)
		{
		ntables = 0;
		tables = NULL;
		checkU_fname = checkF_fname = NULL;
		tables = nullptr;
		}

		/* ---------------------------------------------------------------------- */
		@@ -468,8 +470,6 @@ DihedralTable::~DihedralTable()
		{
		for (int m = 0; m < ntables; m++) free_table(&tables[m]);
		memory->sfree(tables);
		memory->sfree(checkU_fname);
		memory->sfree(checkF_fname);

		if (allocated) {
		memory->destroy(setflag);
		@@ -940,7 +940,7 @@ void DihedralTable::coeff(int narg, char **arg)
		// Optional: allow the user to print out the interpolated spline tables

		if (me == 0) {
		if (checkU_fname && (strlen(checkU_fname) != 0))
		if (!checkU_fname.empty())
		{
		ofstream checkU_file;
		checkU_file.open(checkU_fname, ios::out);
		@@ -953,7 +953,7 @@ void DihedralTable::coeff(int narg, char **arg)
		}
		checkU_file.close();
		}
		if (checkF_fname && (strlen(checkF_fname) != 0))
		if (!checkF_fname.empty())
		{
		ofstream checkF_file;
		checkF_file.open(checkF_fname, ios::out);
		@@ -1084,42 +1084,21 @@ void DihedralTable::free_table(Table *tb)
		/* ----------------------------------------------------------------------
		read table file, only called by proc 0
		------------------------------------------------------------------------- */
		static const int MAXLINE=2048;
		void DihedralTable::read_table(Table tb, char file, char *keyword)
		{
		char line[MAXLINE];
		TableFileReader reader(lmp, file, "dihedral");

		// open file
		char * line = reader.find_section_start(keyword);

		FILE *fp = force->open_potential(file);
		if (fp == NULL) {
		string err_msg = string("Cannot open file ") + string(file);
		error->one(FLERR,err_msg.c_str());
		if (!line) {
		error->one(FLERR,"Did not find keyword in table file");
		}

		// loop until section found with matching keyword

		while (1) {
		if (fgets(line,MAXLINE,fp) == NULL) {
		string err_msg=string("Did not find keyword \"")
		+string(keyword)+string("\" in dihedral table file.");
		error->one(FLERR, err_msg.c_str());
		}
		if (strspn(line," \t\n\r") == strlen(line)) continue; // blank line
		if (line[0] == '#') continue; // comment
		char *word = strtok(line," \t\n\r");
		if (strcmp(word,keyword) == 0) break; // matching keyword
		utils::sfgets(FLERR,line,MAXLINE,fp,file,error); // no match, skip section
		param_extract(tb,line);
		utils::sfgets(FLERR,line,MAXLINE,fp,file,error);
		for (int i = 0; i < tb->ninput; i++)
		utils::sfgets(FLERR,line,MAXLINE,fp,file,error);
		}

		// read args on 2nd line of section
		// allocate table arrays for file values

		utils::sfgets(FLERR,line,MAXLINE,fp,file,error);
		line = reader.next_line();
		param_extract(tb,line);
		memory->create(tb->phifile,tb->ninput,"dihedral:phifile");
		memory->create(tb->efile,tb->ninput,"dihedral:efile");
		@@ -1127,41 +1106,28 @@ void DihedralTable::read_table(Table tb, char file, char *keyword)

		// read a,e,f table values from file

		int itmp;
		for (int i = 0; i < tb->ninput; i++) {
		utils::sfgets(FLERR,line,MAXLINE,fp,file,error);

		// Skip blank lines and delete text following a '#' character
		char *pe = strchr(line, '#');
		if (pe != NULL) *pe = '\0'; //terminate string at '#' character
		char *pc = line;
		while ((pc != '\0') && isspace(pc))
		pc++;
		if (*pc != '\0') { //If line is not a blank line
		stringstream line_ss(line);
		try {
		if (tb->f_unspecified) {
		line_ss >> itmp;
		line_ss >> tb->phifile[i];
		line_ss >> tb->efile[i];
		line = reader.next_line(3);

		ValueTokenizer values(line);
		values.next_int();
		tb->phifile[i] = values.next_double();
		tb->efile[i] = values.next_double();
		} else {
		line_ss >> itmp;
		line_ss >> tb->phifile[i];
		line_ss >> tb->efile[i];
		line_ss >> tb->ffile[i];
		}
		if (! line_ss) {
		stringstream err_msg;
		err_msg << "Read error in table "<< keyword<<", near line "<<i+1<<"\n"
		<< " (Check to make sure the number of columns is correct.)";
		if ((! tb->f_unspecified) && (i==0))
		err_msg << "\n (This sometimes occurs if users forget to specify the \"NOF\" option.)\n";
		error->one(FLERR, err_msg.str().c_str());
		}
		} else //if it is a blank line, then skip it.
		i--;
		} //for (int i = 0; (i < tb->ninput) && fp; i++) {
		line = reader.next_line(4);

		fclose(fp);
		ValueTokenizer values(line);
		values.next_int();
		tb->phifile[i] = values.next_double();
		tb->efile[i] = values.next_double();
		tb->ffile[i] = values.next_double();
		}
		} catch (TokenizerException & e) {
		error->one(FLERR, e.what());
		}
		} //for (int i = 0; (i < tb->ninput) && fp; i++) {
		}

		/* ----------------------------------------------------------------------
		@@ -1353,44 +1319,40 @@ void DihedralTable::param_extract(Table tb, char line)
		tb->f_unspecified = false; //default
		tb->use_degrees = true; //default

		char *word = strtok(line," \t\n\r\f");
		while (word) {
		if (strcmp(word,"N") == 0) {
		word = strtok(NULL," \t\n\r\f");
		tb->ninput = atoi(word);
		try {
		ValueTokenizer values(line);

		while (values.has_next()) {
		std::string word = values.next_string();
		if (word == "N") {
		tb->ninput = values.next_int();
		}
		else if (strcmp(word,"NOF") == 0) {
		else if (word == "NOF") {
		tb->f_unspecified = true;
		}
		else if ((strcmp(word,"DEGREES") == 0) \|\| (strcmp(word,"degrees") == 0)) {
		else if ((word == "DEGREES") \|\| (word == "degrees")) {
		tb->use_degrees = true;
		}
		else if ((strcmp(word,"RADIANS") == 0) \|\| (strcmp(word,"radians") == 0)) {
		else if ((word == "RADIANS") \|\| (word == "radians")) {
		tb->use_degrees = false;
		}
		else if (strcmp(word,"CHECKU") == 0) {
		word = strtok(NULL," \t\n\r\f");
		memory->sfree(checkU_fname);
		memory->create(checkU_fname,strlen(word)+1,"dihedral_table:checkU");
		strcpy(checkU_fname, word);
		else if (word == "CHECKU") {
		checkU_fname = values.next_string();
		}
		else if (strcmp(word,"CHECKF") == 0) {
		word = strtok(NULL," \t\n\r\f");
		memory->sfree(checkF_fname);
		memory->create(checkF_fname,strlen(word)+1,"dihedral_table:checkF");
		strcpy(checkF_fname, word);
		else if (word == "CHECKF") {
		checkF_fname = values.next_string();
		}
		// COMMENTING OUT: equilibrium angles are not supported
		//else if (strcmp(word,"EQ") == 0) {
		// word = strtok(NULL," \t\n\r\f");
		// tb->theta0 = atof(word);
		//else if (word == "EQ") {
		// tb->theta0 = values.next_double();
		//}
		else {
		string err_msg("Invalid keyword in dihedral angle table parameters");
		err_msg += string(" (") + string(word) + string(")");
		string err_msg = fmt::format("Invalid keyword in dihedral angle table parameters ({})", word);
		error->one(FLERR,err_msg.c_str());
		}
		word = strtok(NULL," \t\n\r\f");
		}
		} catch (TokenizerException & e) {
		error->one(FLERR, e.what());
		}

		if (tb->ninput == 0)

src/USER-MISC/dihedral_table.h

+3 −2

Original line number	Diff line number	Diff line
		@@ -24,6 +24,7 @@ DihedralStyle(table,DihedralTable)
		#ifndef LMP_DIHEDRAL_TABLE_H
		#define LMP_DIHEDRAL_TABLE_H
		#include "dihedral.h"
		#include <string>

		namespace LAMMPS_NS {

		@@ -43,8 +44,8 @@ class DihedralTable : public Dihedral {
		protected:
		int tabstyle,tablength;
		// double *phi0; <- equilibrium angles not supported
		char *checkU_fname;
		char *checkF_fname;
		std::string checkU_fname;
		std::string checkF_fname;

		struct Table {
		int ninput;

Admin message