[brahms-md] r49 committed - Added functionalities to generate LAMMPS-compatible input files. Remov...
brah...@googlecode.com
Author: orsimario
Date: Wed Jun 29 10:11:43 2011
Log: Added functionalities to generate LAMMPS-compatible input files.
Removed obsolete references to 'forceNbdStky.c'from Makefiles.
Added 'lammpsScripts' folder to collect LAMMPS-related post-processing
tools.
http://code.google.com/p/brahms-md/source/detail?r=49
Added:
/trunk/lammpsScripts
/trunk/lammpsScripts/lammpstrj2pdb.py
Modified:
/trunk/src/definitions.h
/trunk/src/forceBdMain.c
/trunk/src/main.c
/trunk/src/makefiles/Makefile_gcc-O3
/trunk/src/makefiles/Makefile_icc-fast
/trunk/src/makefiles/Makefile_pgcc-fast
/trunk/src/output.c
=======================================
--- /dev/null
+++ /trunk/lammpsScripts/lammpstrj2pdb.py Wed Jun 29 10:11:43 2011
@@ -0,0 +1,86 @@
+# This script reads the LAMMPS dump file "dump.lammpstrj" and outputs the
+# trajectory file "trajectory.pdb" which can be read in VMD.
+#
+# Dipoles in "dump.lammpstrj" are described by mass center coordinates and
+# orientations. Since VMD does not (currently) allow visualization of
dipoles,
+# this script uses the orientation information to add an extra site to each
+# atom; such an extra site corresponds to the "+" end of the dipole
+#
+# USAGE: python lammpstrj2pdb.py
+
+import sys, string, linecache
+from math import sqrt
+
+print "Processing file dump.lammpstrj ..."
+line = linecache.getline('dump.lammpstrj', 4)
+words = string.split(line)
+nSites = int(words[0])
+print "Number of sites: %d" % nSites
+
+# find the magnitude of the dipole
+line = linecache.getline('dump.lammpstrj', 10)
+words = string.split(line)
+mux = float(words[5])
+muy = float(words[6])
+muz = float(words[7])
+mu = (mux**2 + muy**2 + muz**2)**0.5
+
+print "Dipole magnitude: %.3f" % mu
+
+outFile = open("trajectory.pdb", "w")
+
+inFile = open("dump.lammpstrj", "r")
+lines = inFile.readlines()
+
+print "Setting scaling factor assuming real units..."
+s=2*mu; # scaling factor
+#s=s/30; # rescale for reduced units - comment this line for real units
+nLine=0; # line counter
+
+print "Start scanning coordinates and orientations..."
+for line in lines:
+ nLine=nLine+1
+ words = string.split(line)
+ if len(words) == 2:
+ if words[1]=="TIMESTEP":
+ if nLine!=1:
+ outFile.write('ENDMDL\n')
+ outFile.write('MODEL\n')
+ if len(words) == 8:
+ n=int(words[0])
+ t=int(words[1]) # type
+ x=float(words[2])
+ y=float(words[3])
+ z=float(words[4])
+ mux=float(words[5])
+ muy=float(words[6])
+ muz=float(words[7])
+ xh=x+s*mux
+ yh=y+s*muy
+ zh=z+s*muz
+ if t == 1: # water type
+ outFile.write('ATOM%7d%9s%6d%12.3f%8.3f%8.3f\n'
+ %(n,"OW WAT",n,x,y,z ))
+ nPlus=n+nSites # index of added hydrogen ('+' end of dipole)
+ outFile.write('ATOM%7d%9s%6d%12.3f%8.3f%8.3f\n'
+ %(nPlus,"HW WAT",n,xh,yh,zh ))
+ elif t == 2: # choline type
+ outFile.write('ATOM%7d%9s%6d%12.3f%8.3f%8.3f\n'
+ %(n,"NC3 UNK",n,x,y,z ))
+ elif t == 3: # phosphate type
+ outFile.write('ATOM%7d%9s%6d%12.3f%8.3f%8.3f\n'
+ %(n,"PO4 UNK",n,x,y,z ))
+ elif t == 4 or t == 5: # glycerol or ester type
+ outFile.write('ATOM%7d%9s%6d%12.3f%8.3f%8.3f\n'
+ %(n,"O UNK",n,x,y,z ))
+ nPlus=n+nSites # index of added carbon ('+' end of dipole)
+ outFile.write('ATOM%7d%9s%6d%12.3f%8.3f%8.3f\n'
+ %(nPlus,"C UNK",n,xh,yh,zh ))
+ elif t == 6: # tail type
+ outFile.write('ATOM%7d%9s%6d%12.3f%8.3f%8.3f\n'
+ %(n,"CH2 UNK",n,x,y,z ))
+ elif t == 7: # amine type
+ outFile.write('ATOM%7d%9s%6d%12.3f%8.3f%8.3f\n'
+ %(n,"NH3 UNK",n,x,y,z ))
+
+print "Done scanning - end of script."
=======================================
--- /trunk/src/definitions.h Thu Jan 6 06:46:39 2011
+++ /trunk/src/definitions.h Wed Jun 29 10:11:43 2011
@@ -56,6 +56,7 @@
#define cm_IN_BRAHMS_UNITS ( cm_IN_m * m_IN_BRAHMS_UNITS )
#define D_IN_BRAHMS_UNITS ( DEBYE_IN_C_PER_m * C_IN_BRAHMS_UNITS *
m_IN_BRAHMS_UNITS ) // 1 Debye in brahms units
#define DIPOLE_D ( 1. / D_IN_BRAHMS_UNITS )
+#define D_IN_eA (0.20822678) // 1 D = 0.20822678 electroUnit *
Angstrom
#define s_IN_BRAHMS_UNITS ( 1e12 ) // 1 second in Brahms time units
#define bar_IN_BRAHMS_UNITS ( 1e5 * N_IN_BRAHMS_UNITS / Sqr
(m_IN_BRAHMS_UNITS) ) // 1 bar in Brahms pressure Units
#define atm_IN_BRAHMS_UNITS ( atm_IN_bar * bar_IN_BRAHMS_UNITS ) // 1
bar in Brahms pressure Units
=======================================
--- /trunk/src/forceBdMain.c Thu Jan 6 06:46:39 2011
+++ /trunk/src/forceBdMain.c Wed Jun 29 10:11:43 2011
@@ -30,10 +30,10 @@
/* Global variables: */
real scaleRefSpringLength; // input parameter: referenceSpringLength_ij =
scaleRefSpringLength * (sigma_i+sigma_j)/2
-real refAngle_PhosGlyEst; /* input parameter: reference angle for the
phosphate-glycerol-ester triplet */
real refAngle_CholPhosGly; /* input parameter: reference angle for
phosphocholine headgroup tilt */
-real angleRigidity; /* input parameter: rigidity of the angle potential
for all angles */
+real refAngle_PhosGlyEst; /* input parameter: reference angle for the
phosphate-glycerol-ester triplet */
real refAngle_cisUnsat; /* reference angle for tail sites corresponding to
triplets involving a cis unsaturated double bond */
+real angleRigidity; /* input parameter: rigidity of the angle potential
for all angles */
real orientRigidity; /* input parameter: orientational rigidity for
dipolar lipid sites */
/* CalcAngleBendInteraction -- calculates angle bending interaction,
adapted from [Rapaport, pp.283-285] */
=======================================
--- /trunk/src/main.c Tue Jun 28 10:57:32 2011
+++ /trunk/src/main.c Wed Jun 29 10:11:43 2011
@@ -226,6 +226,7 @@
void UnwrapLipLatCom();
void WriteAreaVolume();
void WriteLammpsData( FILE* );
+void WriteLammpsForcefield( FILE* );
/***********************************************************************
* SetParams -- sets various parameters; adapted from [Rapaport, p.31] *
@@ -281,6 +282,10 @@
SetupJob( fPtr );
fclose( fPtr );
+ fPtr = fopen( "forcefield.elba", "w" );
+ WriteLammpsForcefield( fPtr );
+ fclose( fPtr );
+
fPtr = fopen( "data.initial", "w" );
WriteLammpsData( fPtr );
fclose( fPtr );
=======================================
--- /trunk/src/makefiles/Makefile_gcc-O3 Thu Jan 6 06:46:39 2011
+++ /trunk/src/makefiles/Makefile_gcc-O3 Wed Jun 29 10:11:43 2011
@@ -44,7 +44,6 @@
forceNbdES.c \
forceNbdLJ.c \
forceNbdMain.c \
- forceNbdStky.c \
in_namelist.c \
stepCalc.c \
stepInt.c \
@@ -71,7 +70,6 @@
forceNbdES.o \
forceNbdLJ.o \
forceNbdMain.o \
- forceNbdStky.o \
in_namelist.o \
stepCalc.o \
stepInt.o \
@@ -172,8 +170,6 @@
forceNbdMain.o: forceNbdMain.c dataStructs.h definitions.h
-forceNbdStky.o: forceNbdStky.c dataStructs.h definitions.h
-
rotation.o: rotation.c dataStructs.h definitions.h
startUp.o: startUp.c dataStructs.h definitions.h
=======================================
--- /trunk/src/makefiles/Makefile_icc-fast Thu Jan 6 06:46:39 2011
+++ /trunk/src/makefiles/Makefile_icc-fast Wed Jun 29 10:11:43 2011
@@ -44,7 +44,6 @@
forceNbdES.c \
forceNbdLJ.c \
forceNbdMain.c \
- forceNbdStky.c \
in_namelist.c \
stepCalc.c \
stepInt.c \
@@ -71,7 +70,6 @@
forceNbdES.o \
forceNbdLJ.o \
forceNbdMain.o \
- forceNbdStky.o \
in_namelist.o \
stepCalc.o \
stepInt.o \
@@ -172,8 +170,6 @@
forceNbdMain.o: forceNbdMain.c dataStructs.h definitions.h
-forceNbdStky.o: forceNbdStky.c dataStructs.h definitions.h
-
rotation.o: rotation.c dataStructs.h definitions.h
startUp.o: startUp.c dataStructs.h definitions.h
=======================================
--- /trunk/src/makefiles/Makefile_pgcc-fast Thu Jan 6 06:46:39 2011
+++ /trunk/src/makefiles/Makefile_pgcc-fast Wed Jun 29 10:11:43 2011
@@ -44,7 +44,6 @@
forceNbdES.c \
forceNbdLJ.c \
forceNbdMain.c \
- forceNbdStky.c \
in_namelist.c \
stepCalc.c \
stepInt.c \
@@ -71,7 +70,6 @@
forceNbdES.o \
forceNbdLJ.o \
forceNbdMain.o \
- forceNbdStky.o \
in_namelist.o \
stepCalc.o \
stepInt.o \
@@ -172,8 +170,6 @@
forceNbdMain.o: forceNbdMain.c dataStructs.h definitions.h
-forceNbdStky.o: forceNbdStky.c dataStructs.h definitions.h
-
rotation.o: rotation.c dataStructs.h definitions.h
startUp.o: startUp.c dataStructs.h definitions.h
=======================================
--- /trunk/src/output.c Tue Jun 28 10:57:32 2011
+++ /trunk/src/output.c Wed Jun 29 10:11:43 2011
@@ -19,13 +19,20 @@
extern const Atom *atom;
extern const LipidStruct dopcStruct;
extern const VecRTwo *xyLipidComTrue, xySoluteComTrue;
-extern const real timeNow, sigWat, watDipole, *mass, *charge, *dipole,
*sigma;
+extern const real timeNow, sigWat, watDipole, *mass, *charge, *dipole,
*sigma, **sig, **eps;
+extern const real springRigidity, angleRigidity, refAngle_CholPhosGly,
refAngle_PhosGlyEst, refAngle_cisUnsat;
extern const Prop totEnergy, potEnergy, systemTemperature, lipTemperature,
watTemperature, solTemperature, pressure, dens;
extern const Prop boxVolume, xyArea, lipVolume, lipArea, surfaceTension,
headGroupDipole, headGroupTilt;
extern const Prop ordParTail, ordParTop, ordParMidUnsat, ordParMidSat,
ordParEnd, zConstraintForce;
-extern const int nSites, nWats, nLipids, nSolutes, stepCount, nAtoms,
zConstraint, nTypes;
+extern const int nSites, nWats, nLipids, nDOPCsDSPCs, nDOPEs, nSolutes,
stepCount, nAtoms, zConstraint, nTypes;
extern FILE *fPtrZCF, *fPtrAV;
+// local variables:
+static const int nBondsPerLipid = 14; // assuming 15-site lipid models
+static const int nBondTypes = 6; // 1_chol-phos, 2_phos-gly, 3_gly-est,
4_est-tail, 5_tail-tail, 6_amine_phos
+static const int nAnglesPerLipid = 13; // assuming 15-site lipid models
+static const int nAngleTypes = 4; // 1_chol(amine)-phos-gly,
2_phos-gly-est, 3_PI_sat, 4_cisUnsat
+
/* external function declarations: */
void RotMatToEuler( real*, real*, real*, RMat* );
void RotMatToQuat( Quat*, RMat );
@@ -370,18 +377,23 @@
fprintf( fPtrAV, "%.3f %.2f\n", xyArea.val, boxVolume.val ); /* [nm^2],
[nm^3] */
}
-/****************************************************************************************************************
- * WriteLammpsData -- writes "data." file that can be used by the LAMMPS
program
-
****************************************************************************************************************/
+/*******************************************************************************************************
+ * WriteLammpsData -- writes "data." file to be read into the LAMMPS
program through read_data command *
+
*******************************************************************************************************/
void WriteLammpsData( FILE *fp )
{
- int n;
-
- fprintf(fp, "Brahms-generated input file for LAMMPS\n");
+ int n, bondNumber, angleNumber, atom1, atom2, atom3;
+
+ fprintf(fp, "LAMMPS 'data.' description - created with Brahms
(www.personal.soton.ac.uk/orsi/brahms)\n");
fprintf(fp, "\n");
fprintf(fp, "%9d atoms\n", nSites);
+ fprintf(fp, "%9d bonds\n", nBondsPerLipid*nLipids);
+ fprintf(fp, "%9d angles\n", nAnglesPerLipid*nLipids);
+ fprintf(fp, "\n");
fprintf(fp, "%9d atom types\n", nTypes);
+ fprintf(fp, "%9d bond types\n", nBondTypes);
+ fprintf(fp, "%9d angle types\n", nAngleTypes);
fprintf(fp, "\n");
fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.x*LENGTH_Angstrom/2,
region.x*LENGTH_Angstrom/2, "xlo xhi");
@@ -389,24 +401,213 @@
fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.z*LENGTH_Angstrom/2,
region.z*LENGTH_Angstrom/2, "zlo zhi");
fprintf(fp, "\n");
- fprintf(fp, " Atoms\n");
- fprintf(fp, "\n");
+ fprintf(fp, "Atoms\n\n");
DO_SITE {
// write atom-ID and atom-type:
fprintf(fp, "%7d%5d", n+1, site[n].type+1 );
// write coordinates (x y z):
- fprintf(fp, "%13.5f%13.5f%13.5f", siteCoords[n].x*LENGTH_Angstrom,
siteCoords[n].y*LENGTH_Angstrom, siteCoords[n].z*LENGTH_Angstrom);
+ fprintf(fp, "%11.4f%11.4f%11.4f", siteCoords[n].x*LENGTH_Angstrom,
siteCoords[n].y*LENGTH_Angstrom, siteCoords[n].z*LENGTH_Angstrom);
// write molecule-ID (this is 0 for water sites):
- fprintf(fp, "%6d", site[n].inLipid );
+ fprintf(fp, "%7d", site[n].inLipid );
// write charge (this is 0.00 for water sites):
fprintf(fp, "%7.2f", charge[site[n].type]);
// write dipole component ('dipole' style):
- fprintf(fp, "%13.5f%13.5f%13.5f",
+ fprintf(fp, "% 10.5f% 10.5f% 10.5f",
siteOrientVec[ n ].x * dipole[site[n].type] * DIPOLE_D * D_IN_eA,
siteOrientVec[ n ].y * dipole[site[n].type] * DIPOLE_D * D_IN_eA,
siteOrientVec[ n ].z * dipole[site[n].type] * DIPOLE_D * D_IN_eA );
// write attributes for 'sphere' style - diameter & atom density:
- fprintf(fp, "%7.1f%7.1f", sigma[site[n].type]*LENGTH_Angstrom,
mass[site[n].type] / (4*PI/3 *
Cube((sigma[site[n].type]/2)*LENGTH_Angstrom)));
+ fprintf(fp, "%5.1f%5.1f", sigma[site[n].type]*LENGTH_Angstrom,
mass[site[n].type] / (4*PI/3 *
Cube((sigma[site[n].type]/2)*LENGTH_Angstrom)));
fprintf(fp, "\n");
}
-}
+ fprintf(fp, "\n");
+
+ bondNumber = 0;
+ fprintf(fp, "Bonds\n\n");
+ DO_SITE {
+ switch ( site[ n ].lipidUnit ) {
+ case WATER: /* break - nonlipid site */
+ case TAIL_A5: // break - this case is considered in TAIL_A4
+ case TAIL_B5: // break - this case is considered in TAIL_B4
+ break;
+ case CHOLINE_OR_AMINE: /* choline***phosphate (DOPC) or
amine***phosphate (DOPE) */
+ if (site[n].type == CHOLINE_TYPE) { /* choline***phosphate (DOPC) */
+ atom1 = n + 1;
+ atom2 = atom1 + 1;
+ bondNumber++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondNumber, 1, atom1, atom2 );
+ } else { /* amine***phosphate (DOPE) */
+ atom1 = n + 1;
+ atom2 = atom1 + 1;
+ bondNumber++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondNumber, 6, atom1, atom2 );
+ }
+ break;
+ case PHOSPHATE: // phosphate***glycerol
+ atom1 = n + 1;
+ atom2 = atom1 + 1;
+ bondNumber++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondNumber, 2, atom1, atom2 );
+ break;
+ case GLYCEROL:
+ // GLYCEROL***ESTER_A:
+ atom1 = n + 1;
+ atom2 = atom1 + 1;
+ bondNumber++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondNumber, 3, atom1, atom2 );
+ // GLYCEROL***ESTER_B:
+ atom2 = atom1 + 7;
+ bondNumber++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondNumber, 3, atom1, atom2 );
+ break;
+ case ESTER_A: // ESTER_A***TAIL_A1
+ case ESTER_B: // ESTER_B***TAIL_B1
+ atom1 = n + 1;
+ atom2 = atom1 + 1;
+ bondNumber++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondNumber, 4, atom1, atom2 );
+ break;
+ case TAIL_A1: // TAIL_A1***TAIL_A2
+ case TAIL_B1: // TAIL_B1***TAIL_B2
+ case TAIL_A2: // TAIL_A2***TAIL_A3
+ case TAIL_B2: // TAIL_B2***TAIL_B3
+ case TAIL_A3: // TAIL_A3***TAIL_A4
+ case TAIL_B3: // TAIL_B3***TAIL_B4
+ case TAIL_A4: // do TAIL_A4***TAIL_A5
+ case TAIL_B4: // do TAIL_B4***TAIL_B5
+ atom1 = n + 1;
+ atom2 = atom1 + 1;
+ bondNumber++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondNumber, 5, atom1, atom2 );
+ break;
+ default:
+ printf("Error: default case reached in WriteLammpsData.
Exiting...\n");
+ exit(0);
+ }
+ }
+ fprintf(fp, "\n");
+
+ angleNumber = 0;
+ fprintf(fp, "Angles\n\n");
+ DO_SITE {
+ switch ( site[ n ].lipidUnit ) { /* identify the site "n"
corresponding to the angle "(n-1)--(n)--(n+1)" */
+ case WATER:
+ case CHOLINE_OR_AMINE:
+ case TAIL_A5:
+ case TAIL_B5:
+ break; /* there is no angle potential centered on these sites */
+ case PHOSPHATE: /* choline_or_amine(1)-phosphate(2)-glycerol(3) */
+ angleNumber++;
+ atom1 = n; // choline or amine
+ atom2 = n+1; // phosphate
+ atom3 = n+2; // glycerol
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleNumber, 1, atom1, atom2, atom3
);
+ break;
+ case GLYCEROL:
+ angleNumber++; // phos-gly-estA
+ atom1 = n; // phosphate
+ atom2 = n+1; // glycerol
+ atom3 = n+2; // ester_a
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleNumber, 2, atom1, atom2, atom3
);
+ angleNumber++; // phos-gly-estB
+ atom3 = n+8; // ester_b
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleNumber, 2, atom1, atom2, atom3
);
+ break;
+ case ESTER_A: /* glycerol(3)-ester_A(4)-hydrocarbon_a1(5) */
+ case TAIL_A1: /* ester_A(4)-hydrocarbon_a1(5)-hydrocarbon_a2(6) */
+ case TAIL_A3: /*
hydrocarbon_a2(6)-hydrocarbon_a3(7)-hydrocarbon_a4(8) */
+ case TAIL_A4: /*
hydrocarbon_a3(7)-hydrocarbon_a4(8)-hydrocarbon_a5(9) */
+ case TAIL_B1: /* ester_B(4)-hydrocarbon_b1(5)-hydrocarbon_b2(6) */
+ case TAIL_B3: /*
hydrocarbon_b2(6)-hydrocarbon_b3(7)-hydrocarbon_b4(8) */
+ case TAIL_B4: /*
hydrocarbon_b3(7)-hydrocarbon_b4(8)-hydrocarbon_b5(9) */
+ angleNumber++;
+ atom1 = n;
+ atom2 = n+1;
+ atom3 = n+2;
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleNumber, 3, atom1, atom2, atom3
);
+ break;
+ case ESTER_B: /* glycerol(3)-ester_B(10)-hydrocarbon_b1(11) --
SPECIAL CASE */
+ angleNumber++;
+ atom1 = n-6;
+ atom2 = n+1;
+ atom3 = n+2;
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleNumber, 3, atom1, atom2, atom3
);
+ break;
+ case TAIL_A2: /*
hydrocarbon_a1(5)-hydrocarbon_a2(6)-hydrocarbon_a3(7) */
+ case TAIL_B2: /*
hydrocarbon_b1(5)-hydrocarbon_b2(6)-hydrocarbon_b3(7) */
+ angleNumber++;
+ atom1 = n;
+ atom2 = n+1;
+ atom3 = n+2;
+ // set cisUnsat angle type for DOPC and DOPE - case for DSPC to be
coded!
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleNumber, 4, atom1, atom2, atom3
);
+ break;
+ default:
+ printf("Error in function WriteLammpsData - default case reached.
Exiting...\n"); exit(0);
+ }
+ }
+}
+
+/**********************************************************************************************
+ * WriteLammpsForcefield -- writes "forcefield." file to be included into
LAMMPS input script *
+
**********************************************************************************************/
+void WriteLammpsForcefield( FILE *fp )
+{
+ int i,j,n;
+ real springRigidity_LAMMPS_units = springRigidity * J_IN_cal /
Sqr(nm_IN_Angstrom);
+ real angleRigidity_LAMMPS_units = angleRigidity * J_IN_cal;
+ real radians_in_degrees = ( 180 / PI ); // conversion factor
+
+ // fprintf(fp, "pair_style hybrid lj/cut 12 lj/cut/coul/cut 12
dipole/cut 9\n\n");
+ fprintf(fp, "pair_style dipole/cut 12\n");
+ fprintf(fp, "bond_style harmonic\n");
+ fprintf(fp, "angle_style cosine/squared\n");
+ fprintf(fp, "\n");
+
+ for ( n = 0; n < nTypes; n++ ) { // write masses
+ fprintf(fp, "mass%5d%9.3f\n", n+1, mass[n]);
+ }
+ fprintf(fp, "\n");
+
+ // pair_coeff for 'dipole/cut' style
+ // syntax: pair_coeff iType jType ijEpsilon ijSigma cutoffLJ_optional
cutoffElectro_optional
+ for ( i = 0; i < nTypes; i++ ) {
+ for ( j = 0; j < nTypes; j++ ) {
+ if ( i <= j ) { // LAMMPS restriction to avoid redundancy
+ // fprintf(fp, "pair_coeff%5d%5d\tdipole/cut%7.3f%7.3f\n", i+1,
j+1, eps[i][j]*J_IN_cal, sig[i][j]*nm_IN_Angstrom); // hybrid style
+ fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n", i+1, j+1,
eps[i][j]*J_IN_cal, sig[i][j]*nm_IN_Angstrom);
+ }
+ }
+ }
+ fprintf(fp, "\n");
+
+ // bond_coeff for 'harmonic' style
+ // syntax: bond_coeff bondType rigidityConstant referenceLength
+ // NOTE: in LAMMPS, this potential is defined as K(r-r0)^2, so K_LAMMPS
= k_brahms/2
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 1,
springRigidity_LAMMPS_units/2,
0.9*sig[CHOLINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 2,
springRigidity_LAMMPS_units/2,
0.9*sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 3,
springRigidity_LAMMPS_units/2,
0.9*sig[GLYCEROL_TYPE][ESTER_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 4,
springRigidity_LAMMPS_units/2,
0.9*sig[ESTER_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 5,
springRigidity_LAMMPS_units/2,
0.9*sig[TAIL_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 6,
springRigidity_LAMMPS_units/2,
0.9*sig[PHOSPHATE_TYPE][AMINE_TYPE]*nm_IN_Angstrom);
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 6,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 7,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 8,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 9,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 10,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 11,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 12,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 13,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 14,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ fprintf(fp, "\n");
+
+ // angle_coeff for 'cosine/squared' style
+ // syntax: angle_coeff angleType rigidityConstant referenceAngle
+ // NOTE: in LAMMPS, this potential is defined as
K[cosTheta-cosTheta0]^2, so K_LAMMPS = k_brahms/2
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", 1,
angleRigidity_LAMMPS_units/2, refAngle_CholPhosGly*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", 2,
angleRigidity_LAMMPS_units/2, refAngle_PhosGlyEst*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", 3,
angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", 4,
angleRigidity_LAMMPS_units/2, refAngle_cisUnsat*radians_in_degrees);
+ /* fprintf(fp, "angle_coeff%5d%7.3f%7.3f\n", ,
angleRigidity_LAMMPS_units/2, refAngle_); */
+ fprintf(fp, "\n");
+}