forked from lijiext/lammps
1948 lines
66 KiB
C++
1948 lines
66 KiB
C++
|
/* -----------------------------------------------------------------------
|
||
|
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||
|
|
www.cs.sandia.gov/~sjplimp/lammps.html
|
||
|
|
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
|
||
|
|
|
||
|
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||
|
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||
|
|
certain rights in this software. This software is distributed under
|
||
|
|
the GNU General Public License.
|
||
|
|
|
||
|
|
See the README file in the top-level LAMMPS directory.
|
||
|
|
------------------------------------------------------------------------ */
|
||
|
|
|
||
|
|
// Convert a LAMMPS binary restart file into an ASCII text data file
|
||
|
|
//
|
||
|
|
// Syntax: restart2data restart-file data-file
|
||
|
|
//
|
||
|
|
// this serial code must be compiled on a platform that can read the binary
|
||
|
|
// restart file since binary formats are not compatible across all platforms
|
||
|
|
// restart-file can have a '%' character to indicate a multiproc restart
|
||
|
|
// file as written by LAMMPS
|
||
|
|
|
||
|
|
#include "stdio.h"
|
||
|
|
#include "stdlib.h"
|
||
|
|
#include "string.h"
|
||
|
|
|
||
|
|
#define MIN(a,b) ((a) < (b) ? (a) : (b))
|
||
|
|
#define MAX(a,b) ((a) > (b) ? (a) : (b))
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// Data class to hold problem
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
class Data {
|
||
|
|
public:
|
||
|
|
|
||
|
|
// global settings
|
||
|
|
|
||
|
|
char *version;
|
||
|
|
int ntimestep;
|
||
|
|
int nprocs;
|
||
|
|
char *unit_style;
|
||
|
|
int dimension;
|
||
|
|
int px,py,pz;
|
||
|
|
int newton_pair,newton_bond;
|
||
|
|
int xperiodic,yperiodic,zperiodic;
|
||
|
|
int boundary[3][2];
|
||
|
|
|
||
|
|
char *atom_style;
|
||
|
|
int style_angle,style_atomic,style_bond,style_charge,style_dipole;
|
||
|
|
int style_dpd,style_full,style_granular,style_hybrid,style_molecular;
|
||
|
|
int charge_allow,molecular;
|
||
|
|
int bonds_allow,angles_allow,dihedrals_allow,impropers_allow;
|
||
|
|
|
||
|
|
int natoms,nbonds,nangles,ndihedrals,nimpropers;
|
||
|
|
int ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
|
||
|
|
int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
|
||
|
|
int mass_require,dipole_require;
|
||
|
|
double xlo,xhi,ylo,yhi,zlo,zhi;
|
||
|
|
double special_lj[4],special_coul[4];
|
||
|
|
|
||
|
|
// force fields
|
||
|
|
|
||
|
|
char *pair_style,*bond_style,*angle_style,*dihedral_style,*improper_style;
|
||
|
|
|
||
|
|
double *pair_buck_A,*pair_buck_rho,*pair_buck_C;
|
||
|
|
double *pair_dpd_a0,*pair_dpd_gamma;
|
||
|
|
double *pair_charmm_epsilon,*pair_charmm_sigma;
|
||
|
|
double *pair_charmm_eps14,*pair_charmm_sigma14;
|
||
|
|
double *pair_class2_epsilon,*pair_class2_sigma;
|
||
|
|
double *pair_lj_epsilon,*pair_lj_sigma;
|
||
|
|
double *pair_ljexpand_epsilon,*pair_ljexpand_sigma,*pair_ljexpand_shift;
|
||
|
|
double *pair_morse_d0,*pair_morse_alpha,*pair_morse_r0;
|
||
|
|
double *pair_soft_start,*pair_soft_stop;
|
||
|
|
double *pair_yukawa_A;
|
||
|
|
|
||
|
|
double *bond_class2_r0,*bond_class2_k2,*bond_class2_k3,*bond_class2_k4;
|
||
|
|
double *bond_fene_k,*bond_fene_r0,*bond_fene_epsilon,*bond_fene_sigma;
|
||
|
|
double *bond_feneexpand_k,*bond_feneexpand_r0;
|
||
|
|
double *bond_feneexpand_epsilon,*bond_feneexpand_sigma;
|
||
|
|
double *bond_feneexpand_shift;
|
||
|
|
double *bond_harmonic_k,*bond_harmonic_r0;
|
||
|
|
double *bond_morse_d0,*bond_morse_alpha,*bond_morse_r0;
|
||
|
|
double *bond_nonlinear_epsilon,*bond_nonlinear_r0,*bond_nonlinear_lamda;
|
||
|
|
|
||
|
|
double *angle_charmm_k,*angle_charmm_theta0;
|
||
|
|
double *angle_charmm_k_ub,*angle_charmm_r_ub;
|
||
|
|
double *angle_class2_theta0;
|
||
|
|
double *angle_class2_k2,*angle_class2_k3,*angle_class2_k4;
|
||
|
|
double *angle_class2_bb_k,*angle_class2_bb_r1,*angle_class2_bb_r2;
|
||
|
|
double *angle_class2_ba_k1,*angle_class2_ba_k2;
|
||
|
|
double *angle_class2_ba_r1,*angle_class2_ba_r2;
|
||
|
|
double *angle_cosine_k;
|
||
|
|
double *angle_cosine_squared_k,*angle_cosine_squared_theta0;
|
||
|
|
double *angle_harmonic_k,*angle_harmonic_theta0;
|
||
|
|
|
||
|
|
double *dihedral_charmm_k,*dihedral_charmm_weight;
|
||
|
|
int *dihedral_charmm_multiplicity,*dihedral_charmm_sign;
|
||
|
|
double *dihedral_class2_k1,*dihedral_class2_k2,*dihedral_class2_k3;
|
||
|
|
double *dihedral_class2_phi1,*dihedral_class2_phi2,*dihedral_class2_phi3;
|
||
|
|
double *dihedral_class2_mbt_f1,*dihedral_class2_mbt_f2;
|
||
|
|
double *dihedral_class2_mbt_f3,*dihedral_class2_mbt_r0;
|
||
|
|
double *dihedral_class2_ebt_f1_1,*dihedral_class2_ebt_f2_1;
|
||
|
|
double *dihedral_class2_ebt_f3_1,*dihedral_class2_ebt_r0_1;
|
||
|
|
double *dihedral_class2_ebt_f1_2,*dihedral_class2_ebt_f2_2;
|
||
|
|
double *dihedral_class2_ebt_f3_2,*dihedral_class2_ebt_r0_2;
|
||
|
|
double *dihedral_class2_at_f1_1,*dihedral_class2_at_f2_1;
|
||
|
|
double *dihedral_class2_at_f3_1,*dihedral_class2_at_theta0_1;
|
||
|
|
double *dihedral_class2_at_f1_2,*dihedral_class2_at_f2_2;
|
||
|
|
double *dihedral_class2_at_f3_2,*dihedral_class2_at_theta0_2;
|
||
|
|
double *dihedral_class2_aat_k;
|
||
|
|
double *dihedral_class2_aat_theta0_1,*dihedral_class2_aat_theta0_2;
|
||
|
|
double *dihedral_class2_bb13_k;
|
||
|
|
double *dihedral_class2_bb13_r10,*dihedral_class2_bb13_r30;
|
||
|
|
double *dihedral_harmonic_k;
|
||
|
|
int *dihedral_harmonic_multiplicity,*dihedral_harmonic_sign;
|
||
|
|
double *dihedral_helix_aphi,*dihedral_helix_bphi,*dihedral_helix_cphi;
|
||
|
|
double *dihedral_multi_a1,*dihedral_multi_a2,*dihedral_multi_a3;
|
||
|
|
double *dihedral_multi_a4,*dihedral_multi_a5;
|
||
|
|
|
||
|
|
double *improper_class2_k0,*improper_class2_chi0;
|
||
|
|
double *improper_class2_aa_k1,*improper_class2_aa_k2,*improper_class2_aa_k3;
|
||
|
|
double *improper_class2_aa_theta0_1,*improper_class2_aa_theta0_2;
|
||
|
|
double *improper_class2_aa_theta0_3;
|
||
|
|
double *improper_cvff_k;
|
||
|
|
int *improper_cvff_sign,*improper_cvff_multiplicity;
|
||
|
|
double *improper_harmonic_k,*improper_harmonic_chi;
|
||
|
|
|
||
|
|
// atom quantities
|
||
|
|
|
||
|
|
int iatoms,ibonds,iangles,idihedrals,iimpropers;
|
||
|
|
|
||
|
|
double *mass,*dipole;
|
||
|
|
double *x,*y,*z,*vx,*vy,*vz;
|
||
|
|
double *phix,*phiy,*phiz,*phivx,*phivy,*phivz;
|
||
|
|
int *tag,*type,*mask,*image;
|
||
|
|
int *molecule;
|
||
|
|
double *q,*mux,*muy,*muz,*radius,*density;
|
||
|
|
int *bond_type,*angle_type,*dihedral_type,*improper_type;
|
||
|
|
int *bond_atom1,*bond_atom2;
|
||
|
|
int *angle_atom1,*angle_atom2,*angle_atom3;
|
||
|
|
int *dihedral_atom1,*dihedral_atom2,*dihedral_atom3,*dihedral_atom4;
|
||
|
|
int *improper_atom1,*improper_atom2,*improper_atom3,*improper_atom4;
|
||
|
|
|
||
|
|
// functions
|
||
|
|
|
||
|
|
Data();
|
||
|
|
void stats();
|
||
|
|
void write(FILE *);
|
||
|
|
};
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// function prototypes
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void header(FILE *, Data &);
|
||
|
|
void set_style(char *, Data &);
|
||
|
|
void groups(FILE *);
|
||
|
|
void masses(FILE *, Data &);
|
||
|
|
void dipoles(FILE *, Data &);
|
||
|
|
void force_fields(FILE *, Data &);
|
||
|
|
void modify(FILE *);
|
||
|
|
void pair(FILE *fp, Data &data, char *style, int flag);
|
||
|
|
void bond(FILE *fp, Data &data);
|
||
|
|
void angle(FILE *fp, Data &data);
|
||
|
|
void dihedral(FILE *fp, Data &data);
|
||
|
|
void improper(FILE *fp, Data &data);
|
||
|
|
int atom(double *, Data &data);
|
||
|
|
|
||
|
|
int read_int(FILE *fp);
|
||
|
|
double read_double(FILE *fp);
|
||
|
|
char *read_char(FILE *fp);
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// main program
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
main (int argc, char **argv)
|
||
|
|
{
|
||
|
|
// syntax error check
|
||
|
|
|
||
|
|
if (argc != 3) {
|
||
|
|
printf("Syntax: restart2data restart-file data-file\n");
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
|
||
|
|
// if restart file contains '%', file = filename with % replaced by "base"
|
||
|
|
// else file = single file
|
||
|
|
|
||
|
|
int multiproc;
|
||
|
|
char *file,*ptr;
|
||
|
|
|
||
|
|
if (ptr = strchr(argv[1],'%')) {
|
||
|
|
multiproc = 1;
|
||
|
|
file = new char[strlen(argv[1]) + 16];
|
||
|
|
*ptr = '\0';
|
||
|
|
sprintf(file,"%s%s%s",argv[1],"base",ptr+1);
|
||
|
|
} else {
|
||
|
|
multiproc = 0;
|
||
|
|
file = argv[1];
|
||
|
|
}
|
||
|
|
|
||
|
|
// open single restart file or base file for multiproc case
|
||
|
|
|
||
|
|
printf("Reading restart file ...\n");
|
||
|
|
FILE *fp = fopen(file,"rb");
|
||
|
|
if (fp == NULL) {
|
||
|
|
printf("ERROR: Cannot open restart file %s\n",file);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
|
||
|
|
// read beginning of restart file
|
||
|
|
|
||
|
|
Data data;
|
||
|
|
|
||
|
|
header(fp,data);
|
||
|
|
groups(fp);
|
||
|
|
if (data.mass_require) masses(fp,data);
|
||
|
|
if (data.dipole_require) dipoles(fp,data);
|
||
|
|
force_fields(fp,data);
|
||
|
|
modify(fp);
|
||
|
|
|
||
|
|
// read atoms from single or multiple restart files
|
||
|
|
|
||
|
|
double *buf = NULL;
|
||
|
|
int n,m;
|
||
|
|
int maxbuf = 0;
|
||
|
|
data.iatoms = data.ibonds = data.iangles =
|
||
|
|
data.idihedrals = data.iimpropers = 0;
|
||
|
|
|
||
|
|
for (int iproc = 0; iproc < data.nprocs; iproc++) {
|
||
|
|
if (multiproc) {
|
||
|
|
fclose(fp);
|
||
|
|
sprintf(file,"%s%d%s",argv[1],iproc,ptr+1);
|
||
|
|
fp = fopen(file,"rb");
|
||
|
|
if (fp == NULL) {
|
||
|
|
printf("ERROR: Cannot open restart file %s\n",file);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
n = read_int(fp);
|
||
|
|
|
||
|
|
if (n > maxbuf) {
|
||
|
|
maxbuf = n;
|
||
|
|
delete [] buf;
|
||
|
|
buf = new double[maxbuf];
|
||
|
|
}
|
||
|
|
|
||
|
|
fread(buf,sizeof(double),n,fp);
|
||
|
|
|
||
|
|
m = 0;
|
||
|
|
while (m < n) m += atom(&buf[m],data);
|
||
|
|
}
|
||
|
|
|
||
|
|
fclose(fp);
|
||
|
|
|
||
|
|
// print out stats
|
||
|
|
|
||
|
|
data.stats();
|
||
|
|
|
||
|
|
// write out data file
|
||
|
|
|
||
|
|
printf("Writing data file ...\n");
|
||
|
|
fp = fopen(argv[2],"w");
|
||
|
|
if (fp == NULL) {
|
||
|
|
printf("ERROR: Cannot open data file %s\n",argv[2]);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
|
||
|
|
data.write(fp);
|
||
|
|
fclose(fp);
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// read header of restart file
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void header(FILE *fp, Data &data)
|
||
|
|
{
|
||
|
|
char *version = "1 Oct 2006";
|
||
|
|
|
||
|
|
int flag;
|
||
|
|
flag = read_int(fp);
|
||
|
|
|
||
|
|
while (flag >= 0) {
|
||
|
|
|
||
|
|
if (flag == 0) {
|
||
|
|
data.version = read_char(fp);
|
||
|
|
if (strcmp(version,data.version) != 0) {
|
||
|
|
char *str = "Restart file version does not match restart2data version";
|
||
|
|
printf("WARNING %s\n",str);
|
||
|
|
printf(" restart2data version = %s\n",version);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else if (flag == 1) data.unit_style = read_char(fp);
|
||
|
|
else if (flag == 2) data.ntimestep = read_int(fp);
|
||
|
|
else if (flag == 3) data.dimension = read_int(fp);
|
||
|
|
else if (flag == 4) data.nprocs = read_int(fp);
|
||
|
|
else if (flag == 5) data.px = read_int(fp);
|
||
|
|
else if (flag == 6) data.py = read_int(fp);
|
||
|
|
else if (flag == 7) data.pz = read_int(fp);
|
||
|
|
else if (flag == 8) data.newton_pair = read_int(fp);
|
||
|
|
else if (flag == 9) data.newton_bond = read_int(fp);
|
||
|
|
else if (flag == 10) data.xperiodic = read_int(fp);
|
||
|
|
else if (flag == 11) data.yperiodic = read_int(fp);
|
||
|
|
else if (flag == 12) data.zperiodic = read_int(fp);
|
||
|
|
else if (flag == 13) data.boundary[0][0] = read_int(fp);
|
||
|
|
else if (flag == 14) data.boundary[0][1] = read_int(fp);
|
||
|
|
else if (flag == 15) data.boundary[1][0] = read_int(fp);
|
||
|
|
else if (flag == 16) data.boundary[1][1] = read_int(fp);
|
||
|
|
else if (flag == 17) data.boundary[2][0] = read_int(fp);
|
||
|
|
else if (flag == 18) data.boundary[2][1] = read_int(fp);
|
||
|
|
|
||
|
|
// if atom_style = hybrid, read additional sub-class arguments
|
||
|
|
|
||
|
|
else if (flag == 19) {
|
||
|
|
data.style_angle = data.style_atomic = data.style_bond =
|
||
|
|
data.style_charge = data.style_dipole = data.style_dpd =
|
||
|
|
data.style_full = data.style_granular =
|
||
|
|
data.style_hybrid = data.style_molecular = 0;
|
||
|
|
|
||
|
|
data.atom_style = read_char(fp);
|
||
|
|
set_style(data.atom_style,data);
|
||
|
|
|
||
|
|
if (strcmp(data.atom_style,"hybrid") == 0) {
|
||
|
|
int nwords = read_int(fp);
|
||
|
|
char *substyle;
|
||
|
|
for (int i = 0; i < nwords; i++) {
|
||
|
|
substyle = read_char(fp);
|
||
|
|
set_style(substyle,data);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
else if (flag == 20) data.natoms = static_cast<int> (read_double(fp));
|
||
|
|
else if (flag == 21) data.ntypes = read_int(fp);
|
||
|
|
else if (flag == 22) data.nbonds = read_int(fp);
|
||
|
|
else if (flag == 23) data.nbondtypes = read_int(fp);
|
||
|
|
else if (flag == 24) data.bond_per_atom = read_int(fp);
|
||
|
|
else if (flag == 25) data.nangles = read_int(fp);
|
||
|
|
else if (flag == 26) data.nangletypes = read_int(fp);
|
||
|
|
else if (flag == 27) data.angle_per_atom = read_int(fp);
|
||
|
|
else if (flag == 28) data.ndihedrals = read_int(fp);
|
||
|
|
else if (flag == 29) data.ndihedraltypes = read_int(fp);
|
||
|
|
else if (flag == 30) data.dihedral_per_atom = read_int(fp);
|
||
|
|
else if (flag == 31) data.nimpropers = read_int(fp);
|
||
|
|
else if (flag == 32) data.nimpropertypes = read_int(fp);
|
||
|
|
else if (flag == 33) data.improper_per_atom = read_int(fp);
|
||
|
|
else if (flag == 34) data.xlo = read_double(fp);
|
||
|
|
else if (flag == 35) data.xhi = read_double(fp);
|
||
|
|
else if (flag == 36) data.ylo = read_double(fp);
|
||
|
|
else if (flag == 37) data.yhi = read_double(fp);
|
||
|
|
else if (flag == 38) data.zlo = read_double(fp);
|
||
|
|
else if (flag == 39) data.zhi = read_double(fp);
|
||
|
|
else if (flag == 40) data.special_lj[1] = read_double(fp);
|
||
|
|
else if (flag == 41) data.special_lj[2] = read_double(fp);
|
||
|
|
else if (flag == 42) data.special_lj[3] = read_double(fp);
|
||
|
|
else if (flag == 43) data.special_coul[1] = read_double(fp);
|
||
|
|
else if (flag == 44) data.special_coul[2] = read_double(fp);
|
||
|
|
else if (flag == 45) data.special_coul[3] = read_double(fp);
|
||
|
|
else {
|
||
|
|
printf("ERROR: Invalid flag in header of restart file %d\n",flag);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
|
||
|
|
flag = read_int(fp);
|
||
|
|
}
|
||
|
|
|
||
|
|
// set flags based on atom_style
|
||
|
|
|
||
|
|
data.mass_require = 0;
|
||
|
|
if (data.style_angle || data.style_atomic || data.style_bond ||
|
||
|
|
data.style_charge || data.style_dipole || data.style_dpd ||
|
||
|
|
data.style_full || data.style_hybrid || data.style_molecular)
|
||
|
|
data.mass_require = 1;
|
||
|
|
|
||
|
|
data.charge_allow = 0;
|
||
|
|
if (data.style_charge || data.style_full || data.style_dipole)
|
||
|
|
data.charge_allow = 1;
|
||
|
|
|
||
|
|
data.dipole_require = 0;
|
||
|
|
if (data.style_dipole) data.dipole_require = 1;
|
||
|
|
|
||
|
|
data.molecular = 0;
|
||
|
|
if (data.style_angle || data.style_bond || data.style_full ||
|
||
|
|
data.style_molecular)
|
||
|
|
data.molecular = 1;
|
||
|
|
|
||
|
|
data.bonds_allow = 0;
|
||
|
|
if (data.style_angle || data.style_bond || data.style_full ||
|
||
|
|
data.style_molecular)
|
||
|
|
data.bonds_allow = 1;
|
||
|
|
|
||
|
|
data.angles_allow = 0;
|
||
|
|
if (data.style_angle || data.style_full || data.style_molecular)
|
||
|
|
data.angles_allow = 1;
|
||
|
|
|
||
|
|
data.dihedrals_allow = 0;
|
||
|
|
if (data.style_full || data.style_molecular) data.dihedrals_allow = 1;
|
||
|
|
|
||
|
|
data.impropers_allow = 0;
|
||
|
|
if (data.style_full || data.style_molecular) data.impropers_allow = 1;
|
||
|
|
}
|
||
|
|
|
||
|
|
void set_style(char *name, Data &data)
|
||
|
|
{
|
||
|
|
if (strcmp(name,"angle") == 0) data.style_angle = 1;
|
||
|
|
else if (strcmp(name,"atomic") == 0) data.style_atomic = 1;
|
||
|
|
else if (strcmp(name,"bond") == 0) data.style_bond = 1;
|
||
|
|
else if (strcmp(name,"charge") == 0) data.style_charge = 1;
|
||
|
|
else if (strcmp(name,"dipole") == 0) data.style_dipole = 1;
|
||
|
|
else if (strcmp(name,"dpd") == 0) data.style_dpd = 1;
|
||
|
|
else if (strcmp(name,"full") == 0) data.style_full = 1;
|
||
|
|
else if (strcmp(name,"granular") == 0) data.style_granular = 1;
|
||
|
|
else if (strcmp(name,"hybrid") == 0) data.style_hybrid = 1;
|
||
|
|
else if (strcmp(name,"molecular") == 0) data.style_molecular = 1;
|
||
|
|
else {
|
||
|
|
printf("ERROR: Unknown atom style %s\n",name);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// read group info from restart file, just ignore it
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void groups(FILE *fp)
|
||
|
|
{
|
||
|
|
int ngroup = read_int(fp);
|
||
|
|
|
||
|
|
int n;
|
||
|
|
char *name;
|
||
|
|
|
||
|
|
for (int i = 0; i < ngroup; i++) {
|
||
|
|
name = read_char(fp);
|
||
|
|
delete [] name;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// read masses from restart file
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void masses(FILE *fp, Data &data)
|
||
|
|
{
|
||
|
|
data.mass = new double[data.ntypes+1];
|
||
|
|
fread(&data.mass[1],sizeof(double),data.ntypes,fp);
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// read dipoles from restart file
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void dipoles(FILE *fp, Data &data)
|
||
|
|
{
|
||
|
|
data.dipole = new double[data.ntypes+1];
|
||
|
|
fread(&data.dipole[1],sizeof(double),data.ntypes,fp);
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// read force-field info from restart file
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void force_fields(FILE *fp, Data &data)
|
||
|
|
{
|
||
|
|
int n;
|
||
|
|
|
||
|
|
data.pair_style = read_char(fp);
|
||
|
|
if (data.pair_style) pair(fp,data,data.pair_style,1);
|
||
|
|
|
||
|
|
if (data.bonds_allow) {
|
||
|
|
data.bond_style = read_char(fp);
|
||
|
|
if (data.bond_style) bond(fp,data);
|
||
|
|
}
|
||
|
|
if (data.angles_allow) {
|
||
|
|
data.angle_style = read_char(fp);
|
||
|
|
if (data.angle_style) angle(fp,data);
|
||
|
|
}
|
||
|
|
if (data.dihedrals_allow) {
|
||
|
|
data.dihedral_style = read_char(fp);
|
||
|
|
if (data.dihedral_style) dihedral(fp,data);
|
||
|
|
}
|
||
|
|
if (data.impropers_allow) {
|
||
|
|
data.improper_style = read_char(fp);
|
||
|
|
if (data.improper_style) improper(fp,data);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// read fix info from restart file, just ignore it
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void modify(FILE *fp)
|
||
|
|
{
|
||
|
|
char *buf;
|
||
|
|
int n;
|
||
|
|
|
||
|
|
// nfix = # of fix entries with state
|
||
|
|
|
||
|
|
int nfix = read_int(fp);
|
||
|
|
|
||
|
|
// read each entry with id string, style string, chunk of data
|
||
|
|
|
||
|
|
for (int i = 0; i < nfix; i++) {
|
||
|
|
buf = read_char(fp); delete [] buf;
|
||
|
|
buf = read_char(fp); delete [] buf;
|
||
|
|
buf = read_char(fp); delete [] buf;
|
||
|
|
}
|
||
|
|
|
||
|
|
// nfix = # of fix entries with peratom info
|
||
|
|
|
||
|
|
int nfix_peratom = read_int(fp);
|
||
|
|
|
||
|
|
// read each entry with id string, style string, maxsize of one atom data
|
||
|
|
|
||
|
|
for (int i = 0; i < nfix_peratom; i++) {
|
||
|
|
buf = read_char(fp); delete [] buf;
|
||
|
|
buf = read_char(fp); delete [] buf;
|
||
|
|
n = read_int(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// read atom info from restart file and store in data struct
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
int atom(double *buf, Data &data)
|
||
|
|
{
|
||
|
|
if (data.iatoms == 0) {
|
||
|
|
data.x = new double[data.natoms];
|
||
|
|
data.y = new double[data.natoms];
|
||
|
|
data.z = new double[data.natoms];
|
||
|
|
data.tag = new int[data.natoms];
|
||
|
|
data.type = new int[data.natoms];
|
||
|
|
data.mask = new int[data.natoms];
|
||
|
|
data.image = new int[data.natoms];
|
||
|
|
data.vx = new double[data.natoms];
|
||
|
|
data.vy = new double[data.natoms];
|
||
|
|
data.vz = new double[data.natoms];
|
||
|
|
|
||
|
|
if (data.charge_allow) data.q = new double[data.natoms];
|
||
|
|
|
||
|
|
if (data.style_dipole) {
|
||
|
|
data.mux = new double[data.natoms];
|
||
|
|
data.muy = new double[data.natoms];
|
||
|
|
data.muz = new double[data.natoms];
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.style_granular) {
|
||
|
|
data.phix = new double[data.natoms];
|
||
|
|
data.phiy = new double[data.natoms];
|
||
|
|
data.phiz = new double[data.natoms];
|
||
|
|
data.phivx = new double[data.natoms];
|
||
|
|
data.phivy = new double[data.natoms];
|
||
|
|
data.phivz = new double[data.natoms];
|
||
|
|
data.radius = new double[data.natoms];
|
||
|
|
data.density = new double[data.natoms];
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.molecular) {
|
||
|
|
data.molecule = new int[data.natoms];
|
||
|
|
|
||
|
|
if (data.bonds_allow) {
|
||
|
|
data.bond_type = new int[data.nbonds];
|
||
|
|
data.bond_atom1 = new int[data.nbonds];
|
||
|
|
data.bond_atom2 = new int[data.nbonds];
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.angles_allow) {
|
||
|
|
data.angle_type = new int[data.nangles];
|
||
|
|
data.angle_atom1 = new int[data.nangles];
|
||
|
|
data.angle_atom2 = new int[data.nangles];
|
||
|
|
data.angle_atom3 = new int[data.nangles];
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.dihedrals_allow) {
|
||
|
|
data.dihedral_type = new int[data.ndihedrals];
|
||
|
|
data.dihedral_atom1 = new int[data.ndihedrals];
|
||
|
|
data.dihedral_atom2 = new int[data.ndihedrals];
|
||
|
|
data.dihedral_atom3 = new int[data.ndihedrals];
|
||
|
|
data.dihedral_atom4 = new int[data.ndihedrals];
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.impropers_allow) {
|
||
|
|
data.improper_type = new int[data.nimpropers];
|
||
|
|
data.improper_atom1 = new int[data.nimpropers];
|
||
|
|
data.improper_atom2 = new int[data.nimpropers];
|
||
|
|
data.improper_atom3 = new int[data.nimpropers];
|
||
|
|
data.improper_atom4 = new int[data.nimpropers];
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
int iatoms = data.iatoms;
|
||
|
|
|
||
|
|
int m = 1;
|
||
|
|
|
||
|
|
data.x[iatoms] = buf[m++];
|
||
|
|
data.y[iatoms] = buf[m++];
|
||
|
|
data.z[iatoms] = buf[m++];
|
||
|
|
data.tag[iatoms] = static_cast<int> (buf[m++]);
|
||
|
|
data.type[iatoms] = static_cast<int> (buf[m++]);
|
||
|
|
data.mask[iatoms] = static_cast<int> (buf[m++]);
|
||
|
|
data.image[iatoms] = static_cast<int> (buf[m++]);
|
||
|
|
data.vx[iatoms] = buf[m++];
|
||
|
|
data.vy[iatoms] = buf[m++];
|
||
|
|
data.vz[iatoms] = buf[m++];
|
||
|
|
|
||
|
|
if (data.charge_allow) data.q[iatoms] = buf[m++];
|
||
|
|
|
||
|
|
if (data.style_dipole) {
|
||
|
|
data.mux[iatoms] = buf[m++];
|
||
|
|
data.muy[iatoms] = buf[m++];
|
||
|
|
data.muz[iatoms] = buf[m++];
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.style_granular) {
|
||
|
|
data.phix[iatoms] = buf[m++];
|
||
|
|
data.phiy[iatoms] = buf[m++];
|
||
|
|
data.phiz[iatoms] = buf[m++];
|
||
|
|
data.phivx[iatoms] = buf[m++];
|
||
|
|
data.phivy[iatoms] = buf[m++];
|
||
|
|
data.phivz[iatoms] = buf[m++];
|
||
|
|
data.radius[iatoms] = buf[m++];
|
||
|
|
data.density[iatoms] = buf[m++];
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.molecular) {
|
||
|
|
data.molecule[iatoms] = static_cast<int> (buf[m++]);
|
||
|
|
|
||
|
|
int type,atom1,atom2,atom3,atom4;
|
||
|
|
|
||
|
|
if (data.bonds_allow) {
|
||
|
|
int n = static_cast<int> (buf[m++]);
|
||
|
|
for (int k = 0; k < n; k++) {
|
||
|
|
type = static_cast<int> (buf[m++]);
|
||
|
|
atom1 = static_cast<int> (buf[m++]);
|
||
|
|
if (data.newton_bond || data.tag[iatoms] < atom1) {
|
||
|
|
data.bond_type[data.ibonds] = type;
|
||
|
|
data.bond_atom1[data.ibonds] = data.tag[iatoms];
|
||
|
|
data.bond_atom2[data.ibonds] = atom1;
|
||
|
|
data.ibonds++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.angles_allow) {
|
||
|
|
int n = static_cast<int> (buf[m++]);
|
||
|
|
for (int k = 0; k < n; k++) {
|
||
|
|
type = static_cast<int> (buf[m++]);
|
||
|
|
atom1 = static_cast<int> (buf[m++]);
|
||
|
|
atom2 = static_cast<int> (buf[m++]);
|
||
|
|
atom3 = static_cast<int> (buf[m++]);
|
||
|
|
if (data.newton_bond || data.tag[iatoms] == atom2) {
|
||
|
|
data.angle_type[data.iangles] = type;
|
||
|
|
data.angle_atom1[data.iangles] = atom1;
|
||
|
|
data.angle_atom2[data.iangles] = atom2;
|
||
|
|
data.angle_atom3[data.iangles] = atom3;
|
||
|
|
data.iangles++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.dihedrals_allow) {
|
||
|
|
int n = static_cast<int> (buf[m++]);
|
||
|
|
for (int k = 0; k < n; k++) {
|
||
|
|
type = static_cast<int> (buf[m++]);
|
||
|
|
atom1 = static_cast<int> (buf[m++]);
|
||
|
|
atom2 = static_cast<int> (buf[m++]);
|
||
|
|
atom3 = static_cast<int> (buf[m++]);
|
||
|
|
atom4 = static_cast<int> (buf[m++]);
|
||
|
|
if (data.newton_bond || data.tag[iatoms] == atom2) {
|
||
|
|
data.dihedral_type[data.idihedrals] = type;
|
||
|
|
data.dihedral_atom1[data.idihedrals] = atom1;
|
||
|
|
data.dihedral_atom2[data.idihedrals] = atom2;
|
||
|
|
data.dihedral_atom3[data.idihedrals] = atom3;
|
||
|
|
data.dihedral_atom4[data.idihedrals] = atom4;
|
||
|
|
data.idihedrals++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (data.impropers_allow) {
|
||
|
|
int n = static_cast<int> (buf[m++]);
|
||
|
|
for (int k = 0; k < n; k++) {
|
||
|
|
type = static_cast<int> (buf[m++]);
|
||
|
|
atom1 = static_cast<int> (buf[m++]);
|
||
|
|
atom2 = static_cast<int> (buf[m++]);
|
||
|
|
atom3 = static_cast<int> (buf[m++]);
|
||
|
|
atom4 = static_cast<int> (buf[m++]);
|
||
|
|
if (data.newton_bond || data.tag[iatoms] == atom2) {
|
||
|
|
data.improper_type[data.iimpropers] = type;
|
||
|
|
data.improper_atom1[data.iimpropers] = atom1;
|
||
|
|
data.improper_atom2[data.iimpropers] = atom2;
|
||
|
|
data.improper_atom3[data.iimpropers] = atom3;
|
||
|
|
data.improper_atom4[data.iimpropers] = atom4;
|
||
|
|
data.iimpropers++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
data.iatoms++;
|
||
|
|
m = static_cast<int> (buf[0]);
|
||
|
|
return m;
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// pair coeffs
|
||
|
|
// one section for each pair style
|
||
|
|
// flag = 1, read all coeff info and allocation arrays
|
||
|
|
// flag = 0, just read global settings (when called recursively by hybrid)
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void pair(FILE *fp, Data &data, char *style, int flag)
|
||
|
|
{
|
||
|
|
int i,j,m;
|
||
|
|
int itmp;
|
||
|
|
double rtmp;
|
||
|
|
|
||
|
|
if (strcmp(style,"none") == 0) {
|
||
|
|
|
||
|
|
} else if ((strcmp(style,"buck") == 0) ||
|
||
|
|
(strcmp(style,"buck/coul/cut") == 0) ||
|
||
|
|
(strcmp(style,"buck/coul/long") == 0)) {
|
||
|
|
|
||
|
|
if (strcmp(style,"buck") == 0) {
|
||
|
|
m = 0;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"buck/coul/cut") == 0) {
|
||
|
|
m = 1;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
double cut_lj_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"buck/coul/long") == 0) {
|
||
|
|
m = 0;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
double cut_lj_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (!flag) return;
|
||
|
|
|
||
|
|
data.pair_buck_A = new double[data.ntypes+1];
|
||
|
|
data.pair_buck_rho = new double[data.ntypes+1];
|
||
|
|
data.pair_buck_C = new double[data.ntypes+1];
|
||
|
|
|
||
|
|
for (i = 1; i <= data.ntypes; i++)
|
||
|
|
for (j = i; j <= data.ntypes; j++) {
|
||
|
|
itmp = read_int(fp);
|
||
|
|
if (i == j && itmp == 0) {
|
||
|
|
printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
if (itmp) {
|
||
|
|
if (i == j) {
|
||
|
|
data.pair_buck_A[i] = read_double(fp);
|
||
|
|
data.pair_buck_rho[i] = read_double(fp);
|
||
|
|
data.pair_buck_C[i] = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
if (m) double cut_coul = read_double(fp);
|
||
|
|
} else {
|
||
|
|
double buck_A = read_double(fp);
|
||
|
|
double buck_rho = read_double(fp);
|
||
|
|
double buck_C = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
if (m) double cut_coul = read_double(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
} else if (strcmp(style,"dpd") == 0) {
|
||
|
|
|
||
|
|
double temperature = read_double(fp);
|
||
|
|
double cut_global = read_double(fp);
|
||
|
|
int seed = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
|
||
|
|
if (!flag) return;
|
||
|
|
|
||
|
|
data.pair_dpd_a0 = new double[data.ntypes+1];
|
||
|
|
data.pair_dpd_gamma = new double[data.ntypes+1];
|
||
|
|
|
||
|
|
for (i = 1; i <= data.ntypes; i++)
|
||
|
|
for (j = i; j <= data.ntypes; j++) {
|
||
|
|
itmp = read_int(fp);
|
||
|
|
if (i == j && itmp == 0) {
|
||
|
|
printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
if (itmp) {
|
||
|
|
data.pair_dpd_a0[i] = read_double(fp);
|
||
|
|
data.pair_dpd_gamma[i] = read_double(fp);
|
||
|
|
double cut = read_double(fp);
|
||
|
|
} else {
|
||
|
|
double dpd_da0 = read_double(fp);
|
||
|
|
double dpd_gamma = read_double(fp);
|
||
|
|
double cut = read_double(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
} else if (strcmp(style,"eam") == 0) {
|
||
|
|
} else if (strcmp(style,"eam/alloy") == 0) {
|
||
|
|
} else if (strcmp(style,"eam/fs") == 0) {
|
||
|
|
|
||
|
|
} else if ((strcmp(style,"gran/history") == 0) ||
|
||
|
|
(strcmp(style,"gran/no_history") == 0) ||
|
||
|
|
(strcmp(style,"gran/hertzian") == 0)) {
|
||
|
|
|
||
|
|
double xkk = read_double(fp);
|
||
|
|
double gamman = read_double(fp);
|
||
|
|
double xmu = read_double(fp);
|
||
|
|
int dampflag = read_int(fp);
|
||
|
|
|
||
|
|
} else if ((strcmp(style,"lj/charmm/coul/charmm") == 0) ||
|
||
|
|
(strcmp(style,"lj/charmm/coul/charmm/implicit") == 0) ||
|
||
|
|
(strcmp(style,"lj/charmm/coul/long") == 0)) {
|
||
|
|
|
||
|
|
if (strcmp(style,"lj/charmm/coul/charmm") == 0) {
|
||
|
|
double cut_lj_inner = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
double cut_coul_inner = read_double(fp);
|
||
|
|
double cut_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"lj/charmm/coul/charmm/implicit") == 0) {
|
||
|
|
double cut_lj_inner = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
double cut_coul_inner = read_double(fp);
|
||
|
|
double cut_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"lj/charmm/coul/long") == 0) {
|
||
|
|
double cut_lj_inner = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
double cut_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (!flag) return;
|
||
|
|
|
||
|
|
data.pair_charmm_epsilon = new double[data.ntypes+1];
|
||
|
|
data.pair_charmm_sigma = new double[data.ntypes+1];
|
||
|
|
data.pair_charmm_eps14 = new double[data.ntypes+1];
|
||
|
|
data.pair_charmm_sigma14 = new double[data.ntypes+1];
|
||
|
|
|
||
|
|
for (i = 1; i <= data.ntypes; i++)
|
||
|
|
for (j = i; j <= data.ntypes; j++) {
|
||
|
|
itmp = read_int(fp);
|
||
|
|
if (i == j && itmp == 0) {
|
||
|
|
printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
if (itmp) {
|
||
|
|
if (i == j) {
|
||
|
|
data.pair_charmm_epsilon[i] = read_double(fp);
|
||
|
|
data.pair_charmm_sigma[i] = read_double(fp);
|
||
|
|
data.pair_charmm_eps14[i] = read_double(fp);
|
||
|
|
data.pair_charmm_sigma14[i] = read_double(fp);
|
||
|
|
} else {
|
||
|
|
double charmm_epsilon = read_double(fp);
|
||
|
|
double charmm_sigma = read_double(fp);
|
||
|
|
double charmm_eps14 = read_double(fp);
|
||
|
|
double charmm_sigma14 = read_double(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
} else if ((strcmp(style,"lj/class2") == 0) ||
|
||
|
|
(strcmp(style,"lj/class2/coul/cut") == 0) ||
|
||
|
|
(strcmp(style,"lj/class2/coul/long") == 0)) {
|
||
|
|
|
||
|
|
if (strcmp(style,"lj/class2") == 0) {
|
||
|
|
m = 0;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"lj/class2/coul/cut") == 0) {
|
||
|
|
m = 1;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
double cut_lj_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"lj/class2/coul/long") == 0) {
|
||
|
|
m = 0;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
double cut_lj_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (!flag) return;
|
||
|
|
|
||
|
|
data.pair_class2_epsilon = new double[data.ntypes+1];
|
||
|
|
data.pair_class2_sigma = new double[data.ntypes+1];
|
||
|
|
|
||
|
|
for (i = 1; i <= data.ntypes; i++)
|
||
|
|
for (j = i; j <= data.ntypes; j++) {
|
||
|
|
itmp = read_int(fp);
|
||
|
|
if (i == j && itmp == 0) {
|
||
|
|
printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
if (itmp) {
|
||
|
|
if (i == j) {
|
||
|
|
data.pair_class2_epsilon[i] = read_double(fp);
|
||
|
|
data.pair_class2_sigma[i] = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
if (m) double cut_coul = read_double(fp);
|
||
|
|
} else {
|
||
|
|
double class2_epsilon = read_double(fp);
|
||
|
|
double class2_sigma = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
if (m) double cut_coul = read_double(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
} else if ((strcmp(style,"lj/cut") == 0) ||
|
||
|
|
(strcmp(style,"lj/cut/coul/cut") == 0) ||
|
||
|
|
(strcmp(style,"lj/cut/coul/debye") == 0) ||
|
||
|
|
(strcmp(style,"lj/cut/coul/long") == 0) ||
|
||
|
|
(strcmp(style,"lj/cut/coul/long/tip4p") == 0)) {
|
||
|
|
|
||
|
|
if (strcmp(style,"lj/cut") == 0) {
|
||
|
|
m = 0;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"lj/cut/coul/cut") == 0) {
|
||
|
|
m = 1;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
double cut_lj_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"lj/cut/coul/debye") == 0) {
|
||
|
|
m = 1;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
double cut_lj_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"lj/cut/coul/long") == 0) {
|
||
|
|
m = 0;
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
double cut_lj_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
} else if (strcmp(style,"lj/cut/coul/long/tip4p") == 0) {
|
||
|
|
m = 0;
|
||
|
|
int typeO = read_int(fp);
|
||
|
|
int typeH = read_int(fp);
|
||
|
|
int typeB = read_int(fp);
|
||
|
|
int typeA = read_int(fp);
|
||
|
|
double qdist = read_double(fp);
|
||
|
|
double cut_lj_global = read_double(fp);
|
||
|
|
double cut_lj_coul = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (!flag) return;
|
||
|
|
|
||
|
|
data.pair_lj_epsilon = new double[data.ntypes+1];
|
||
|
|
data.pair_lj_sigma = new double[data.ntypes+1];
|
||
|
|
|
||
|
|
for (i = 1; i <= data.ntypes; i++)
|
||
|
|
for (j = i; j <= data.ntypes; j++) {
|
||
|
|
itmp = read_int(fp);
|
||
|
|
if (i == j && itmp == 0) {
|
||
|
|
printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
if (itmp) {
|
||
|
|
if (i == j) {
|
||
|
|
data.pair_lj_epsilon[i] = read_double(fp);
|
||
|
|
data.pair_lj_sigma[i] = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
if (m) double cut_coul = read_double(fp);
|
||
|
|
} else {
|
||
|
|
double lj_epsilon = read_double(fp);
|
||
|
|
double lj_sigma = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
if (m) double cut_coul = read_double(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
} else if (strcmp(style,"lj/expand") == 0) {
|
||
|
|
|
||
|
|
double cut_global = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
|
||
|
|
if (!flag) return;
|
||
|
|
|
||
|
|
data.pair_ljexpand_epsilon = new double[data.ntypes+1];
|
||
|
|
data.pair_ljexpand_sigma = new double[data.ntypes+1];
|
||
|
|
data.pair_ljexpand_shift = new double[data.ntypes+1];
|
||
|
|
|
||
|
|
for (i = 1; i <= data.ntypes; i++)
|
||
|
|
for (j = i; j <= data.ntypes; j++) {
|
||
|
|
itmp = read_int(fp);
|
||
|
|
if (i == j && itmp == 0) {
|
||
|
|
printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
if (itmp) {
|
||
|
|
data.pair_ljexpand_epsilon[i] = read_double(fp);
|
||
|
|
data.pair_ljexpand_sigma[i] = read_double(fp);
|
||
|
|
data.pair_ljexpand_shift[i] = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
} else {
|
||
|
|
double ljexpand_epsilon = read_double(fp);
|
||
|
|
double ljexpand_sigma = read_double(fp);
|
||
|
|
double ljexpand_shift = read_double(fp);
|
||
|
|
double cut_lj = read_double(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
} else if (strcmp(style,"morse") == 0) {
|
||
|
|
|
||
|
|
double cut_global = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
|
||
|
|
if (!flag) return;
|
||
|
|
|
||
|
|
data.pair_morse_d0 = new double[data.ntypes+1];
|
||
|
|
data.pair_morse_alpha = new double[data.ntypes+1];
|
||
|
|
data.pair_morse_r0 = new double[data.ntypes+1];
|
||
|
|
|
||
|
|
for (i = 1; i <= data.ntypes; i++)
|
||
|
|
for (j = i; j <= data.ntypes; j++) {
|
||
|
|
itmp = read_int(fp);
|
||
|
|
if (i == j && itmp == 0) {
|
||
|
|
printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
if (itmp) {
|
||
|
|
data.pair_morse_d0[i] = read_double(fp);
|
||
|
|
data.pair_morse_alpha[i] = read_double(fp);
|
||
|
|
data.pair_morse_r0[i] = read_double(fp);
|
||
|
|
double cut = read_double(fp);
|
||
|
|
} else {
|
||
|
|
double morse_d0 = read_double(fp);
|
||
|
|
double morse_alpha = read_double(fp);
|
||
|
|
double morse_r0 = read_double(fp);
|
||
|
|
double cut = read_double(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
} else if (strcmp(style,"soft") == 0) {
|
||
|
|
|
||
|
|
double cut_global = read_double(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
|
||
|
|
if (!flag) return;
|
||
|
|
|
||
|
|
data.pair_soft_start = new double[data.ntypes+1];
|
||
|
|
data.pair_soft_stop = new double[data.ntypes+1];
|
||
|
|
|
||
|
|
for (i = 1; i <= data.ntypes; i++)
|
||
|
|
for (j = i; j <= data.ntypes; j++) {
|
||
|
|
itmp = read_int(fp);
|
||
|
|
if (i == j && itmp == 0) {
|
||
|
|
printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
if (itmp) {
|
||
|
|
data.pair_soft_start[i] = read_double(fp);
|
||
|
|
data.pair_soft_stop[i] = read_double(fp);
|
||
|
|
double cut = read_double(fp);
|
||
|
|
} else {
|
||
|
|
double soft_start = read_double(fp);
|
||
|
|
double soft_stop = read_double(fp);
|
||
|
|
double cut = read_double(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
} else if (strcmp(style,"table") == 0) {
|
||
|
|
|
||
|
|
int tabstyle = read_int(fp);
|
||
|
|
int n = read_int(fp);
|
||
|
|
|
||
|
|
} else if (strcmp(style,"yukawa") == 0) {
|
||
|
|
|
||
|
|
double kappa = read_double(fp);
|
||
|
|
double cut_global = read_double(fp);
|
||
|
|
int offset_flag = read_int(fp);
|
||
|
|
int mix_flag = read_int(fp);
|
||
|
|
|
||
|
|
if (!flag) return;
|
||
|
|
|
||
|
|
data.pair_yukawa_A = new double[data.ntypes+1];
|
||
|
|
|
||
|
|
for (i = 1; i <= data.ntypes; i++)
|
||
|
|
for (j = i; j <= data.ntypes; j++) {
|
||
|
|
itmp = read_int(fp);
|
||
|
|
if (i == j && itmp == 0) {
|
||
|
|
printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
if (itmp) {
|
||
|
|
data.pair_yukawa_A[i] = read_double(fp);
|
||
|
|
double cut = read_double(fp);
|
||
|
|
} else {
|
||
|
|
double yukawa_A = read_double(fp);
|
||
|
|
double cut = read_double(fp);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
} else if (strcmp(style,"hybrid") == 0) {
|
||
|
|
|
||
|
|
// for each substyle of hybrid,
|
||
|
|
// read its settings by calling pair() recursively with flag = 0
|
||
|
|
// so that coeff array allocation is skipped
|
||
|
|
|
||
|
|
int nstyles = read_int(fp);
|
||
|
|
for (int i = 0; i < nstyles; i++) {
|
||
|
|
char *substyle = read_char(fp);
|
||
|
|
pair(fp,data,substyle,0);
|
||
|
|
}
|
||
|
|
|
||
|
|
} else {
|
||
|
|
printf("ERROR: Unknown pair style %s\n",style);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// bond coeffs
|
||
|
|
// one section for each bond style
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void bond(FILE *fp, Data &data)
|
||
|
|
{
|
||
|
|
if (strcmp(data.bond_style,"none") == 0) {
|
||
|
|
|
||
|
|
} else if (strcmp(data.bond_style,"class2") == 0) {
|
||
|
|
|
||
|
|
data.bond_class2_r0 = new double[data.nbondtypes+1];
|
||
|
|
data.bond_class2_k2 = new double[data.nbondtypes+1];
|
||
|
|
data.bond_class2_k3 = new double[data.nbondtypes+1];
|
||
|
|
data.bond_class2_k4 = new double[data.nbondtypes+1];
|
||
|
|
fread(&data.bond_class2_r0[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_class2_k2[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_class2_k3[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_class2_k4[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.bond_style,"fene") == 0) {
|
||
|
|
|
||
|
|
data.bond_fene_k = new double[data.nbondtypes+1];
|
||
|
|
data.bond_fene_r0 = new double[data.nbondtypes+1];
|
||
|
|
data.bond_fene_epsilon = new double[data.nbondtypes+1];
|
||
|
|
data.bond_fene_sigma = new double[data.nbondtypes+1];
|
||
|
|
fread(&data.bond_fene_k[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_fene_r0[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_fene_epsilon[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_fene_sigma[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.bond_style,"fene/expand") == 0) {
|
||
|
|
|
||
|
|
data.bond_feneexpand_k = new double[data.nbondtypes+1];
|
||
|
|
data.bond_feneexpand_r0 = new double[data.nbondtypes+1];
|
||
|
|
data.bond_feneexpand_epsilon = new double[data.nbondtypes+1];
|
||
|
|
data.bond_feneexpand_sigma = new double[data.nbondtypes+1];
|
||
|
|
data.bond_feneexpand_shift = new double[data.nbondtypes+1];
|
||
|
|
fread(&data.bond_feneexpand_k[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_feneexpand_r0[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_feneexpand_epsilon[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_feneexpand_sigma[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_feneexpand_shift[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.bond_style,"harmonic") == 0) {
|
||
|
|
|
||
|
|
data.bond_harmonic_k = new double[data.nbondtypes+1];
|
||
|
|
data.bond_harmonic_r0 = new double[data.nbondtypes+1];
|
||
|
|
fread(&data.bond_harmonic_k[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_harmonic_r0[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.bond_style,"morse") == 0) {
|
||
|
|
|
||
|
|
data.bond_morse_d0 = new double[data.nbondtypes+1];
|
||
|
|
data.bond_morse_alpha = new double[data.nbondtypes+1];
|
||
|
|
data.bond_morse_r0 = new double[data.nbondtypes+1];
|
||
|
|
fread(&data.bond_morse_d0[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_morse_alpha[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_morse_r0[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.bond_style,"nonlinear") == 0) {
|
||
|
|
|
||
|
|
data.bond_nonlinear_epsilon = new double[data.nbondtypes+1];
|
||
|
|
data.bond_nonlinear_r0 = new double[data.nbondtypes+1];
|
||
|
|
data.bond_nonlinear_lamda = new double[data.nbondtypes+1];
|
||
|
|
fread(&data.bond_nonlinear_epsilon[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_nonlinear_r0[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
fread(&data.bond_nonlinear_lamda[1],sizeof(double),data.nbondtypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.bond_style,"hybrid") == 0) {
|
||
|
|
|
||
|
|
int nstyles = read_int(fp);
|
||
|
|
for (int i = 0; i < nstyles; i++)
|
||
|
|
char *substyle = read_char(fp);
|
||
|
|
|
||
|
|
} else {
|
||
|
|
printf("ERROR: Unknown bond style %s\n",data.bond_style);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// angle coeffs
|
||
|
|
// one section for each angle style
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void angle(FILE *fp, Data &data)
|
||
|
|
{
|
||
|
|
if (strcmp(data.angle_style,"none") == 0) {
|
||
|
|
|
||
|
|
} else if (strcmp(data.angle_style,"charmm") == 0) {
|
||
|
|
|
||
|
|
data.angle_charmm_k = new double[data.nangletypes+1];
|
||
|
|
data.angle_charmm_theta0 = new double[data.nangletypes+1];
|
||
|
|
data.angle_charmm_k_ub = new double[data.nangletypes+1];
|
||
|
|
data.angle_charmm_r_ub = new double[data.nangletypes+1];
|
||
|
|
fread(&data.angle_charmm_k[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_charmm_theta0[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_charmm_k_ub[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_charmm_r_ub[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.angle_style,"class2") == 0) {
|
||
|
|
|
||
|
|
data.angle_class2_theta0 = new double[data.nangletypes+1];
|
||
|
|
data.angle_class2_k2 = new double[data.nangletypes+1];
|
||
|
|
data.angle_class2_k3 = new double[data.nangletypes+1];
|
||
|
|
data.angle_class2_k4 = new double[data.nangletypes+1];
|
||
|
|
|
||
|
|
data.angle_class2_bb_k = new double[data.nangletypes+1];
|
||
|
|
data.angle_class2_bb_r1 = new double[data.nangletypes+1];
|
||
|
|
data.angle_class2_bb_r2 = new double[data.nangletypes+1];
|
||
|
|
|
||
|
|
data.angle_class2_ba_k1 = new double[data.nangletypes+1];
|
||
|
|
data.angle_class2_ba_k2 = new double[data.nangletypes+1];
|
||
|
|
data.angle_class2_ba_r1 = new double[data.nangletypes+1];
|
||
|
|
data.angle_class2_ba_r2 = new double[data.nangletypes+1];
|
||
|
|
|
||
|
|
fread(&data.angle_class2_theta0[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_class2_k2[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_class2_k3[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_class2_k4[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
|
||
|
|
fread(&data.angle_class2_bb_k[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_class2_bb_r1[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_class2_bb_r2[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
|
||
|
|
fread(&data.angle_class2_ba_k1[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_class2_ba_k2[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_class2_ba_r1[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_class2_ba_r2[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.angle_style,"cosine") == 0) {
|
||
|
|
|
||
|
|
data.angle_cosine_k = new double[data.nangletypes+1];
|
||
|
|
fread(&data.angle_cosine_k[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.angle_style,"cosine/squared") == 0) {
|
||
|
|
|
||
|
|
data.angle_cosine_squared_k = new double[data.nangletypes+1];
|
||
|
|
data.angle_cosine_squared_theta0 = new double[data.nangletypes+1];
|
||
|
|
fread(&data.angle_cosine_squared_k[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_cosine_squared_theta0[1],
|
||
|
|
sizeof(double),data.nangletypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.angle_style,"harmonic") == 0) {
|
||
|
|
|
||
|
|
data.angle_harmonic_k = new double[data.nangletypes+1];
|
||
|
|
data.angle_harmonic_theta0 = new double[data.nangletypes+1];
|
||
|
|
fread(&data.angle_harmonic_k[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
fread(&data.angle_harmonic_theta0[1],sizeof(double),data.nangletypes,fp);
|
||
|
|
|
||
|
|
} else {
|
||
|
|
printf("ERROR: Unknown angle style %s\n",data.angle_style);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// dihedral coeffs
|
||
|
|
// one section for each dihedral style
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void dihedral(FILE *fp, Data &data)
|
||
|
|
{
|
||
|
|
if (strcmp(data.dihedral_style,"none") == 0) {
|
||
|
|
|
||
|
|
} else if (strcmp(data.dihedral_style,"charmm") == 0) {
|
||
|
|
|
||
|
|
data.dihedral_charmm_k = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_charmm_multiplicity = new int[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_charmm_sign = new int[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_charmm_weight = new double[data.ndihedraltypes+1];
|
||
|
|
fread(&data.dihedral_charmm_k[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_charmm_multiplicity[1],sizeof(int),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_charmm_sign[1],sizeof(int),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_charmm_weight[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.dihedral_style,"class2") == 0) {
|
||
|
|
|
||
|
|
data.dihedral_class2_k1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_k2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_k3 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_phi1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_phi2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_phi3 = new double[data.ndihedraltypes+1];
|
||
|
|
|
||
|
|
data.dihedral_class2_mbt_f1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_mbt_f2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_mbt_f3 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_mbt_r0 = new double[data.ndihedraltypes+1];
|
||
|
|
|
||
|
|
data.dihedral_class2_ebt_f1_1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_ebt_f2_1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_ebt_f3_1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_ebt_r0_1 = new double[data.ndihedraltypes+1];
|
||
|
|
|
||
|
|
data.dihedral_class2_ebt_f1_2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_ebt_f2_2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_ebt_f3_2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_ebt_r0_2 = new double[data.ndihedraltypes+1];
|
||
|
|
|
||
|
|
data.dihedral_class2_at_f1_1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_at_f2_1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_at_f3_1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_at_theta0_1 = new double[data.ndihedraltypes+1];
|
||
|
|
|
||
|
|
data.dihedral_class2_at_f1_2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_at_f2_2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_at_f3_2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_at_theta0_2 = new double[data.ndihedraltypes+1];
|
||
|
|
|
||
|
|
data.dihedral_class2_aat_k = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_aat_theta0_1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_aat_theta0_2 = new double[data.ndihedraltypes+1];
|
||
|
|
|
||
|
|
data.dihedral_class2_bb13_k = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_bb13_r10 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_class2_bb13_r30 = new double[data.ndihedraltypes+1];
|
||
|
|
|
||
|
|
fread(&data.dihedral_class2_k1[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_k2[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_k3[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_phi1[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_phi2[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_phi3[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
fread(&data.dihedral_class2_mbt_f1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_mbt_f2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_mbt_f3[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_mbt_r0[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
fread(&data.dihedral_class2_ebt_f1_1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_ebt_f2_1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_ebt_f3_1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_ebt_r0_1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
fread(&data.dihedral_class2_ebt_f1_2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_ebt_f2_2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_ebt_f3_2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_ebt_r0_2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
fread(&data.dihedral_class2_at_f1_1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_at_f2_1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_at_f3_1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_at_theta0_1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
fread(&data.dihedral_class2_at_f1_2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_at_f2_2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_at_f3_2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_at_theta0_2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
fread(&data.dihedral_class2_aat_k[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_aat_theta0_1[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_aat_theta0_2[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
fread(&data.dihedral_class2_bb13_k[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_bb13_r10[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_class2_bb13_r30[1],sizeof(double),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.dihedral_style,"harmonic") == 0) {
|
||
|
|
|
||
|
|
data.dihedral_harmonic_k = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_harmonic_multiplicity = new int[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_harmonic_sign = new int[data.ndihedraltypes+1];
|
||
|
|
fread(&data.dihedral_harmonic_k[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_harmonic_multiplicity[1],sizeof(int),
|
||
|
|
data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_harmonic_sign[1],sizeof(int),data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.dihedral_style,"helix") == 0) {
|
||
|
|
|
||
|
|
data.dihedral_helix_aphi = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_helix_bphi = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_helix_cphi = new double[data.ndihedraltypes+1];
|
||
|
|
fread(&data.dihedral_helix_aphi[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_helix_bphi[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_helix_cphi[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.dihedral_style,"multi/harmonic") == 0) {
|
||
|
|
|
||
|
|
data.dihedral_multi_a1 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_multi_a2 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_multi_a3 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_multi_a4 = new double[data.ndihedraltypes+1];
|
||
|
|
data.dihedral_multi_a5 = new double[data.ndihedraltypes+1];
|
||
|
|
fread(&data.dihedral_multi_a1[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_multi_a2[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_multi_a3[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_multi_a4[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
fread(&data.dihedral_multi_a5[1],sizeof(double),data.ndihedraltypes,fp);
|
||
|
|
|
||
|
|
} else {
|
||
|
|
printf("ERROR: Unknown dihedral style %s\n",data.dihedral_style);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// improper coeffs
|
||
|
|
// one section for each improper style
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void improper(FILE *fp, Data &data)
|
||
|
|
{
|
||
|
|
if (strcmp(data.improper_style,"none") == 0) {
|
||
|
|
|
||
|
|
} else if (strcmp(data.improper_style,"class2") == 0) {
|
||
|
|
|
||
|
|
data.improper_class2_k0 = new double[data.nimpropertypes+1];
|
||
|
|
data.improper_class2_chi0 = new double[data.nimpropertypes+1];
|
||
|
|
|
||
|
|
data.improper_class2_aa_k1 = new double[data.nimpropertypes+1];
|
||
|
|
data.improper_class2_aa_k2 = new double[data.nimpropertypes+1];
|
||
|
|
data.improper_class2_aa_k3 = new double[data.nimpropertypes+1];
|
||
|
|
data.improper_class2_aa_theta0_1 = new double[data.nimpropertypes+1];
|
||
|
|
data.improper_class2_aa_theta0_2 = new double[data.nimpropertypes+1];
|
||
|
|
data.improper_class2_aa_theta0_3 = new double[data.nimpropertypes+1];
|
||
|
|
|
||
|
|
fread(&data.improper_class2_k0[1],sizeof(double),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
fread(&data.improper_class2_chi0[1],sizeof(double),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
|
||
|
|
fread(&data.improper_class2_aa_k1[1],sizeof(double),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
fread(&data.improper_class2_aa_k2[1],sizeof(double),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
fread(&data.improper_class2_aa_k3[1],sizeof(double),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
fread(&data.improper_class2_aa_theta0_1[1],sizeof(double),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
fread(&data.improper_class2_aa_theta0_2[1],sizeof(double),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
fread(&data.improper_class2_aa_theta0_3[1],sizeof(double),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.improper_style,"cvff") == 0) {
|
||
|
|
|
||
|
|
data.improper_cvff_k = new double[data.nimpropertypes+1];
|
||
|
|
data.improper_cvff_sign = new int[data.nimpropertypes+1];
|
||
|
|
data.improper_cvff_multiplicity = new int[data.nimpropertypes+1];
|
||
|
|
fread(&data.improper_cvff_k[1],sizeof(double),data.nimpropertypes,fp);
|
||
|
|
fread(&data.improper_cvff_sign[1],sizeof(int),data.nimpropertypes,fp);
|
||
|
|
fread(&data.improper_cvff_multiplicity[1],sizeof(int),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
|
||
|
|
} else if (strcmp(data.improper_style,"harmonic") == 0) {
|
||
|
|
|
||
|
|
data.improper_harmonic_k = new double[data.nimpropertypes+1];
|
||
|
|
data.improper_harmonic_chi = new double[data.nimpropertypes+1];
|
||
|
|
fread(&data.improper_harmonic_k[1],sizeof(double),data.nimpropertypes,fp);
|
||
|
|
fread(&data.improper_harmonic_chi[1],sizeof(double),
|
||
|
|
data.nimpropertypes,fp);
|
||
|
|
|
||
|
|
} else {
|
||
|
|
printf("ERROR: Unknown improper style %s\n",data.improper_style);
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// initialize ptrs in Data
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
Data::Data()
|
||
|
|
{
|
||
|
|
pair_style = bond_style = angle_style =
|
||
|
|
dihedral_style = improper_style = NULL;
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// print out stats on problem
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void Data::stats()
|
||
|
|
{
|
||
|
|
printf(" Restart file version = %s\n",version);
|
||
|
|
printf(" Ntimestep = %d\n",ntimestep);
|
||
|
|
printf(" Nprocs = %d\n",nprocs);
|
||
|
|
printf(" Natoms = %d\n",natoms);
|
||
|
|
printf(" Nbonds = %d\n",nbonds);
|
||
|
|
printf(" Nangles = %d\n",nangles);
|
||
|
|
printf(" Ndihedrals = %d\n",ndihedrals);
|
||
|
|
printf(" Nimpropers = %d\n",nimpropers);
|
||
|
|
printf(" Unit style = %s\n",unit_style);
|
||
|
|
printf(" Atom style = %s\n",atom_style);
|
||
|
|
printf(" Pair style = %s\n",pair_style);
|
||
|
|
printf(" Bond style = %s\n",bond_style);
|
||
|
|
printf(" Angle style = %s\n",angle_style);
|
||
|
|
printf(" Dihedral style = %s\n",dihedral_style);
|
||
|
|
printf(" Improper style = %s\n",improper_style);
|
||
|
|
printf(" Xlo/xhi = %g %g\n",xlo,xhi);
|
||
|
|
printf(" Ylo/yhi = %g %g\n",ylo,yhi);
|
||
|
|
printf(" Zlo/zhi = %g %g\n",zlo,zhi);
|
||
|
|
printf(" Periodicity = %d %d %d\n",xperiodic,yperiodic,zperiodic);
|
||
|
|
printf(" Boundary = %d %d, %d %d, %d %d\n",boundary[0][0],boundary[0][1],
|
||
|
|
boundary[1][0],boundary[1][1],boundary[2][0],boundary[2][1]);
|
||
|
|
}
|
||
|
|
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
// write the data file
|
||
|
|
// ---------------------------------------------------------------------
|
||
|
|
|
||
|
|
void Data::write(FILE *fp)
|
||
|
|
{
|
||
|
|
fprintf(fp,"LAMMPS data file from restart file: timestep = %d, procs = %d\n",
|
||
|
|
ntimestep,nprocs);
|
||
|
|
|
||
|
|
fprintf(fp,"\n");
|
||
|
|
|
||
|
|
fprintf(fp,"%d atoms\n",natoms);
|
||
|
|
if (nbonds) fprintf(fp,"%d bonds\n",nbonds);
|
||
|
|
if (nangles) fprintf(fp,"%d angles\n",nangles);
|
||
|
|
if (ndihedrals) fprintf(fp,"%d dihedrals\n",ndihedrals);
|
||
|
|
if (nimpropers) fprintf(fp,"%d impropers\n",nimpropers);
|
||
|
|
|
||
|
|
fprintf(fp,"\n");
|
||
|
|
|
||
|
|
fprintf(fp,"%d atom types\n",ntypes);
|
||
|
|
if (nbondtypes) fprintf(fp,"%d bond types\n",nbondtypes);
|
||
|
|
if (nangletypes) fprintf(fp,"%d angle types\n",nangletypes);
|
||
|
|
if (ndihedraltypes) fprintf(fp,"%d dihedral types\n",ndihedraltypes);
|
||
|
|
if (nimpropertypes) fprintf(fp,"%d improper types\n",nimpropertypes);
|
||
|
|
|
||
|
|
fprintf(fp,"\n");
|
||
|
|
|
||
|
|
fprintf(fp,"%g %g xlo xhi\n",xlo,xhi);
|
||
|
|
fprintf(fp,"%g %g ylo yhi\n",ylo,yhi);
|
||
|
|
fprintf(fp,"%g %g zlo zhi\n",zlo,zhi);
|
||
|
|
|
||
|
|
if (mass_require) {
|
||
|
|
fprintf(fp,"\nMasses\n\n");
|
||
|
|
for (int i = 1; i <= ntypes; i++) fprintf(fp,"%d %g\n",i,mass[i]);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (dipole_require) {
|
||
|
|
fprintf(fp,"\nDipoles\n\n");
|
||
|
|
for (int i = 1; i <= ntypes; i++) fprintf(fp,"%d %g\n",i,dipole[i]);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (pair_style) {
|
||
|
|
if ((strcmp(pair_style,"none") != 0) &&
|
||
|
|
(strcmp(pair_style,"eam") != 0) &&
|
||
|
|
(strcmp(pair_style,"gran/history") != 0) &&
|
||
|
|
(strcmp(pair_style,"gran/no_history") != 0) &&
|
||
|
|
(strcmp(pair_style,"gran/hertzian") != 0) &&
|
||
|
|
(strcmp(pair_style,"table") != 0) &&
|
||
|
|
(strcmp(pair_style,"hybrid") != 0))
|
||
|
|
fprintf(fp,"\nPair Coeffs\n\n");
|
||
|
|
|
||
|
|
if ((strcmp(pair_style,"buck") == 0) ||
|
||
|
|
(strcmp(pair_style,"buck/coul/cut") == 0) ||
|
||
|
|
(strcmp(pair_style,"buck/coul/long") == 0)) {
|
||
|
|
for (int i = 1; i <= ntypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g\n",i,
|
||
|
|
pair_buck_A[i],pair_buck_rho[i],pair_buck_C[i]);
|
||
|
|
|
||
|
|
} else if ((strcmp(pair_style,"lj/charmm/coul/charmm") == 0) ||
|
||
|
|
(strcmp(pair_style,"lj/charmm/coul/long") == 0) ||
|
||
|
|
(strcmp(pair_style,"lj/charmm/coul/charmm/implicit") == 0)) {
|
||
|
|
for (int i = 1; i <= ntypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g\n",i,
|
||
|
|
pair_charmm_epsilon[i],pair_charmm_sigma[i],
|
||
|
|
pair_charmm_eps14[i],pair_charmm_sigma14[i]);
|
||
|
|
|
||
|
|
} else if ((strcmp(pair_style,"lj/class2") == 0) ||
|
||
|
|
(strcmp(pair_style,"lj/class2/coul/cut") == 0) ||
|
||
|
|
(strcmp(pair_style,"lj/class2/coul/long") == 0)) {
|
||
|
|
for (int i = 1; i <= ntypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g\n",i,
|
||
|
|
pair_class2_epsilon[i],pair_class2_sigma[i]);
|
||
|
|
|
||
|
|
} else if ((strcmp(pair_style,"lj/cut") == 0) ||
|
||
|
|
(strcmp(pair_style,"lj/cut/coul/cut") == 0) ||
|
||
|
|
(strcmp(pair_style,"lj/cut/coul/debye") == 0) ||
|
||
|
|
(strcmp(pair_style,"lj/cut/coul/long") == 0) ||
|
||
|
|
(strcmp(pair_style,"lj/cut/coul/long/tip4p") == 0)) {
|
||
|
|
for (int i = 1; i <= ntypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g\n",i,
|
||
|
|
pair_lj_epsilon[i],pair_lj_sigma[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(pair_style,"lj/expand") == 0) {
|
||
|
|
for (int i = 1; i <= ntypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g\n",i,
|
||
|
|
pair_ljexpand_epsilon[i],pair_ljexpand_sigma[i],
|
||
|
|
pair_ljexpand_shift[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(pair_style,"morse") == 0) {
|
||
|
|
for (int i = 1; i <= ntypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g\n",i,
|
||
|
|
pair_morse_d0[i],pair_morse_alpha[i],pair_morse_r0[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(pair_style,"soft") == 0) {
|
||
|
|
for (int i = 1; i <= ntypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g\n",i,
|
||
|
|
pair_soft_start[i],pair_soft_stop[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(pair_style,"yukawa") == 0) {
|
||
|
|
for (int i = 1; i <= ntypes; i++)
|
||
|
|
fprintf(fp,"%d %g\n",i,
|
||
|
|
pair_yukawa_A[i]);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (bond_style) {
|
||
|
|
if ((strcmp(bond_style,"none") != 0) &&
|
||
|
|
(strcmp(bond_style,"hybrid") != 0))
|
||
|
|
fprintf(fp,"\nBond Coeffs\n\n");
|
||
|
|
|
||
|
|
if (strcmp(bond_style,"class2") == 0) {
|
||
|
|
for (int i = 1; i <= nbondtypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g\n",i,
|
||
|
|
bond_class2_r0[i],bond_class2_k2[i],
|
||
|
|
bond_class2_k3[i],bond_class2_k4[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(bond_style,"fene") == 0) {
|
||
|
|
for (int i = 1; i <= nbondtypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g\n",i,
|
||
|
|
bond_fene_k[i],bond_fene_r0[i],
|
||
|
|
bond_fene_epsilon[i],bond_fene_sigma[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(bond_style,"fene/expand") == 0) {
|
||
|
|
for (int i = 1; i <= nbondtypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g %g\n",i,
|
||
|
|
bond_feneexpand_k[i],bond_feneexpand_r0[i],
|
||
|
|
bond_feneexpand_epsilon[i],bond_feneexpand_sigma[i],
|
||
|
|
bond_feneexpand_shift[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(bond_style,"harmonic") == 0) {
|
||
|
|
for (int i = 1; i <= nbondtypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g\n",i,
|
||
|
|
bond_harmonic_k[i],bond_harmonic_r0[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(bond_style,"morse") == 0) {
|
||
|
|
for (int i = 1; i <= nbondtypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g\n",i,
|
||
|
|
bond_morse_d0[i],bond_morse_alpha[i],bond_morse_r0[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(bond_style,"nonlinear") == 0) {
|
||
|
|
for (int i = 1; i <= nbondtypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g\n",i,
|
||
|
|
bond_nonlinear_epsilon[i],bond_nonlinear_r0[i],
|
||
|
|
bond_nonlinear_lamda[i]);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (angle_style) {
|
||
|
|
double PI = 3.1415926; // convert back to degrees
|
||
|
|
|
||
|
|
if (strcmp(angle_style,"none") != 0)
|
||
|
|
fprintf(fp,"\nAngle Coeffs\n\n");
|
||
|
|
|
||
|
|
if (strcmp(angle_style,"charmm") == 0) {
|
||
|
|
for (int i = 1; i <= nangletypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g\n",i,
|
||
|
|
angle_charmm_k[i],angle_charmm_theta0[i]/PI*180.0,
|
||
|
|
angle_charmm_k_ub[i],angle_charmm_r_ub[i]);
|
||
|
|
|
||
|
|
} if (strcmp(angle_style,"class2") == 0) {
|
||
|
|
for (int i = 1; i <= nangletypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g\n",i,
|
||
|
|
angle_class2_theta0[i]/PI*180.0,angle_class2_k2[i],
|
||
|
|
angle_class2_k3[i],angle_class2_k4[i]);
|
||
|
|
|
||
|
|
fprintf(fp,"\nBondBond Coeffs\n\n");
|
||
|
|
for (int i = 1; i <= nangletypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g\n",i,
|
||
|
|
angle_class2_bb_k[i],
|
||
|
|
angle_class2_bb_r1[i],angle_class2_bb_r2[i]);
|
||
|
|
|
||
|
|
fprintf(fp,"\nBondAngle Coeffs\n\n");
|
||
|
|
for (int i = 1; i <= nangletypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g\n",i,
|
||
|
|
angle_class2_ba_k1[i],angle_class2_ba_k2[i],
|
||
|
|
angle_class2_ba_r1[i],angle_class2_ba_r2[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(angle_style,"cosine") == 0) {
|
||
|
|
for (int i = 1; i <= nangletypes; i++)
|
||
|
|
fprintf(fp,"%d %g\n",i,
|
||
|
|
angle_cosine_k[i]);
|
||
|
|
|
||
|
|
} if (strcmp(angle_style,"harmonic") == 0) {
|
||
|
|
for (int i = 1; i <= nangletypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g\n",i,
|
||
|
|
angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (dihedral_style) {
|
||
|
|
double PI = 3.1415926; // convert back to degrees
|
||
|
|
|
||
|
|
if (strcmp(dihedral_style,"none") != 0)
|
||
|
|
fprintf(fp,"\nDihedral Coeffs\n\n");
|
||
|
|
|
||
|
|
if (strcmp(dihedral_style,"charmm") == 0) {
|
||
|
|
for (int i = 1; i <= ndihedraltypes; i++)
|
||
|
|
fprintf(fp,"%d %g %d %d %g\n",i,
|
||
|
|
dihedral_charmm_k[i],dihedral_charmm_multiplicity[i],
|
||
|
|
dihedral_charmm_sign[i],dihedral_charmm_weight[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(dihedral_style,"class2") == 0) {
|
||
|
|
for (int i = 1; i <= ndihedraltypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g %g %g\n",i,
|
||
|
|
dihedral_class2_k1[i],dihedral_class2_k2[i],
|
||
|
|
dihedral_class2_k3[i],
|
||
|
|
dihedral_class2_phi1[i]/PI*180.0,
|
||
|
|
dihedral_class2_phi2[i]/PI*180.0,
|
||
|
|
dihedral_class2_phi3[i]/PI*180.0);
|
||
|
|
|
||
|
|
fprintf(fp,"\nMiddleBondTorsion Coeffs\n\n");
|
||
|
|
for (int i = 1; i <= ndihedraltypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g\n",i,
|
||
|
|
dihedral_class2_mbt_f1[i],dihedral_class2_mbt_f2[i],
|
||
|
|
dihedral_class2_mbt_f3[i],dihedral_class2_mbt_r0[i]);
|
||
|
|
|
||
|
|
fprintf(fp,"\nEndBondTorsion Coeffs\n\n");
|
||
|
|
for (int i = 1; i <= ndihedraltypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g %g %g %g %g\n",i,
|
||
|
|
dihedral_class2_ebt_f1_1[i],dihedral_class2_ebt_f2_1[i],
|
||
|
|
dihedral_class2_ebt_f3_1[i],dihedral_class2_ebt_r0_1[i],
|
||
|
|
dihedral_class2_ebt_f1_2[i],dihedral_class2_ebt_f2_2[i],
|
||
|
|
dihedral_class2_ebt_f3_2[i],dihedral_class2_ebt_r0_2[i]);
|
||
|
|
|
||
|
|
fprintf(fp,"\nAngleTorsion Coeffs\n\n");
|
||
|
|
for (int i = 1; i <= ndihedraltypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g %g %g %g %g\n",i,
|
||
|
|
dihedral_class2_at_f1_1[i],dihedral_class2_at_f2_1[i],
|
||
|
|
dihedral_class2_at_f3_1[i],
|
||
|
|
dihedral_class2_at_theta0_1[i]/PI*180.0,
|
||
|
|
dihedral_class2_at_f1_2[i],dihedral_class2_at_f2_2[i],
|
||
|
|
dihedral_class2_at_f3_2[i],
|
||
|
|
dihedral_class2_at_theta0_2[i]/PI*180.0);
|
||
|
|
|
||
|
|
fprintf(fp,"\nAngleAngleTorsion Coeffs\n\n");
|
||
|
|
for (int i = 1; i <= ndihedraltypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g\n",i,
|
||
|
|
dihedral_class2_aat_k[i],
|
||
|
|
dihedral_class2_aat_theta0_1[i]/PI*180.0,
|
||
|
|
dihedral_class2_aat_theta0_2[i]/PI*180.0);
|
||
|
|
|
||
|
|
fprintf(fp,"\nBondBond13 Coeffs\n\n");
|
||
|
|
for (int i = 1; i <= ndihedraltypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g\n",i,
|
||
|
|
dihedral_class2_bb13_k[i],
|
||
|
|
dihedral_class2_bb13_r10[i],dihedral_class2_bb13_r30[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(dihedral_style,"harmonic") == 0) {
|
||
|
|
for (int i = 1; i <= ndihedraltypes; i++)
|
||
|
|
fprintf(fp,"%d %g %d %d\n",i,
|
||
|
|
dihedral_harmonic_k[i],dihedral_harmonic_multiplicity[i],
|
||
|
|
dihedral_harmonic_sign[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(dihedral_style,"multi/harmonic") == 0) {
|
||
|
|
for (int i = 1; i <= ndihedraltypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g %g\n",i,
|
||
|
|
dihedral_multi_a1[i],dihedral_multi_a2[i],
|
||
|
|
dihedral_multi_a3[i],dihedral_multi_a4[i],
|
||
|
|
dihedral_multi_a5[i]);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (improper_style) {
|
||
|
|
double PI = 3.1415926; // convert back to degrees
|
||
|
|
|
||
|
|
if (strcmp(improper_style,"none") != 0)
|
||
|
|
fprintf(fp,"\nImproper Coeffs\n\n");
|
||
|
|
|
||
|
|
if (strcmp(improper_style,"class2") == 0) {
|
||
|
|
for (int i = 1; i <= nimpropertypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g\n",i,
|
||
|
|
improper_class2_k0[i],improper_class2_chi0[i]/PI*180.0);
|
||
|
|
|
||
|
|
fprintf(fp,"\nAngleAngle Coeffs\n\n");
|
||
|
|
for (int i = 1; i <= nimpropertypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g %g %g %g %g\n",i,
|
||
|
|
improper_class2_aa_k1[i],improper_class2_aa_k2[i],
|
||
|
|
improper_class2_aa_k3[i],
|
||
|
|
improper_class2_aa_theta0_1[i]/PI*180.0,
|
||
|
|
improper_class2_aa_theta0_2[i]/PI*180.0,
|
||
|
|
improper_class2_aa_theta0_3[i]/PI*180.0);
|
||
|
|
|
||
|
|
} else if (strcmp(improper_style,"cvff") == 0) {
|
||
|
|
for (int i = 1; i <= nimpropertypes; i++)
|
||
|
|
fprintf(fp,"%d %g %d %d\n",i,
|
||
|
|
improper_cvff_k[i],improper_cvff_sign[i],
|
||
|
|
improper_cvff_multiplicity[i]);
|
||
|
|
|
||
|
|
} else if (strcmp(improper_style,"harmonic") == 0) {
|
||
|
|
for (int i = 1; i <= nimpropertypes; i++)
|
||
|
|
fprintf(fp,"%d %g %g\n",i,
|
||
|
|
improper_harmonic_k[i],improper_harmonic_chi[i]/PI*180.0);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (natoms) {
|
||
|
|
fprintf(fp,"\nAtoms\n\n");
|
||
|
|
|
||
|
|
int ix,iy,iz;
|
||
|
|
for (int i = 0; i < natoms; i++) {
|
||
|
|
ix = image[i]%1000 - 500;
|
||
|
|
iy = (image[i]/1000)%1000 - 500;
|
||
|
|
iz = image[i]/1000000 - 500;
|
||
|
|
|
||
|
|
fprintf(fp,"%d",tag[i]);
|
||
|
|
if (molecular) fprintf(fp," %d",molecule[i]);
|
||
|
|
fprintf(fp," %d",type[i]);
|
||
|
|
if (charge_allow) fprintf(fp," %g",q[i]);
|
||
|
|
if (style_granular) fprintf(fp," %g %g",2.0*radius[i],density[i]);
|
||
|
|
fprintf(fp," %g %g %g",x[i],y[i],z[i]);
|
||
|
|
if (style_dipole) fprintf(fp," %g %g %g",mux[i],muy[i],muz[i]);
|
||
|
|
fprintf(fp," %d %d %d\n",ix,iy,iz);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (natoms) {
|
||
|
|
fprintf(fp,"\nVelocities\n\n");
|
||
|
|
for (int i = 0; i < natoms; i++)
|
||
|
|
if (strcmp(atom_style,"granular") == 0)
|
||
|
|
fprintf(fp,"%d %g %g %g %g %g %g\n",tag[i],vx[i],vy[i],vz[i],
|
||
|
|
phivx[i],phivy[i],phivz[i]);
|
||
|
|
else
|
||
|
|
fprintf(fp,"%d %g %g %g\n",tag[i],vx[i],vy[i],vz[i]);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (nbonds) {
|
||
|
|
fprintf(fp,"\nBonds\n\n");
|
||
|
|
for (int i = 0; i < nbonds; i++)
|
||
|
|
fprintf(fp,"%d %d %d %d\n",
|
||
|
|
i+1,bond_type[i],bond_atom1[i],bond_atom2[i]);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (nangles) {
|
||
|
|
fprintf(fp,"\nAngles\n\n");
|
||
|
|
for (int i = 0; i < nangles; i++)
|
||
|
|
fprintf(fp,"%d %d %d %d %d\n",
|
||
|
|
i+1,angle_type[i],angle_atom1[i],angle_atom2[i],angle_atom3[i]);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (ndihedrals) {
|
||
|
|
fprintf(fp,"\nDihedrals\n\n");
|
||
|
|
for (int i = 0; i < ndihedrals; i++)
|
||
|
|
fprintf(fp,"%d %d %d %d %d %d\n",
|
||
|
|
i+1,dihedral_type[i],dihedral_atom1[i],dihedral_atom2[i],
|
||
|
|
dihedral_atom3[i],dihedral_atom4[i]);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (nimpropers) {
|
||
|
|
fprintf(fp,"\nImpropers\n\n");
|
||
|
|
for (int i = 0; i < nimpropers; i++)
|
||
|
|
fprintf(fp,"%d %d %d %d %d %d\n",
|
||
|
|
i+1,improper_type[i],improper_atom1[i],improper_atom2[i],
|
||
|
|
improper_atom3[i],improper_atom4[i]);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
/* read an int from restart file */
|
||
|
|
|
||
|
|
int read_int(FILE *fp)
|
||
|
|
{
|
||
|
|
int value;
|
||
|
|
fread(&value,sizeof(int),1,fp);
|
||
|
|
return value;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* read a double from restart file */
|
||
|
|
|
||
|
|
double read_double(FILE *fp)
|
||
|
|
{
|
||
|
|
double value;
|
||
|
|
fread(&value,sizeof(double),1,fp);
|
||
|
|
return value;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* read a char str from restart file */
|
||
|
|
|
||
|
|
char *read_char(FILE *fp)
|
||
|
|
{
|
||
|
|
int n;
|
||
|
|
fread(&n,sizeof(int),1,fp);
|
||
|
|
if (n == 0) return NULL;
|
||
|
|
char *value = new char[n];
|
||
|
|
fread(value,sizeof(char),n,fp);
|
||
|
|
return value;
|
||
|
|
}
|