[brahms-md] r50 committed - More work on Brahms/LAMMPS interfacing

0 views

Skip to first unread message

brah...@googlecode.com

unread,

Jul 6, 2011, 1:05:05 PM7/6/11

to brahms-d...@googlegroups.com

Revision: 50
Author: orsimario
Date: Wed Jul 6 10:04:36 2011
Log: More work on Brahms/LAMMPS interfacing
http://code.google.com/p/brahms-md/source/detail?r=50

Added:
/trunk/src/lammps.c
Modified:
/trunk/lammpsScripts/lammpstrj2pdb.py
/trunk/src/definitions.h
/trunk/src/makefiles/Makefile_gcc-O3
/trunk/src/makefiles/Makefile_icc-fast
/trunk/src/output.c

=======================================
--- /dev/null
+++ /trunk/src/lammps.c Wed Jul 6 10:04:36 2011
@@ -0,0 +1,381 @@
+/* Copyright (C) 2011 Mario Orsi
+ This file is part of Brahms.
+ Brahms is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.
+ Brahms is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
+ You should have received a copy of the GNU General Public License along
with Brahms. If not, see <http://www.gnu.org/licenses/>. */
+
+/***********************************************************
+ * lammps.c -- functions for interfacing Brahms and LAMMPS *
+ ***********************************************************/
+
+#include "dataStructs.h"
+
+extern const VecR region, *siteCoords, *siteOrientVec;
+extern const Site *site;
+extern const real timeNow, *mass, *charge, *dipole, *sigma, **sig, **eps;
+extern const real springRigidity, angleRigidity, refAngle_CholPhosGly,
refAngle_PhosGlyEst, refAngle_cisUnsat;
+extern const int nSites, nLipids, nDOPCsDSPCs, nDOPEs, nTypes, nWaters;
+
+// local variables:
+static const int nAtomsPerLipid = 21; // assuming 15-site lipid models
with 2 added atoms for each gly and est site
+static const int nBondsPerLipid = 23; // assuming 15-site lipid models: 14
(main) + 6 (tip-tail charges) + 3 (orientation-restraining)
+static const int nBondTypes = 9; // 1_chol-phos, 2_phos-gly, 3_gly-est,
4_est-tail, 5_tail-tail, 6_amine_phos, 7_glyPlus_phos, 8_estPlus_tail,
9_atomicCharges
+static const int nAnglesPerLipid = 16; // assuming 15-site lipid models:
13 (main) + 3 (atomic charges)
+static const int nAngleTypes = 5; // 1_chol(amine)-phos-gly,
2_phos-gly-est, 3_PI_sat, 4_cisUnsat, 5_atomicCharges
+
+/*******************************************************************************************************
+ * WriteLammpsData -- writes "data." file to be read into the LAMMPS
program through read_data command *
+
*******************************************************************************************************/
+void WriteLammpsData( FILE *fp )
+{
+ int n, id, atomCount, bondCount, angleCount, atom1, atom2, atom3;
+ int estCom, estDipTail, estDipTip, chain;
+ int nAtoms = nWaters + nAtomsPerLipid*nLipids; // #atoms for LAMMPS
topology
+ int nBonds = nBondsPerLipid*nLipids;
+ int nAngles = nAnglesPerLipid*nLipids;
+ real disp, atomCharge;
+ VecR tailCoords, tipCoords;
+
+ fprintf(fp, "LAMMPS 'data.' description - created with Brahms
(www.personal.soton.ac.uk/orsi/brahms)\n");
+ fprintf(fp, "\n");
+
+ fprintf(fp, "%9d atoms\n", nAtoms);
+ fprintf(fp, "%9d bonds\n", nBonds);
+ fprintf(fp, "%9d angles\n", nAngles);
+ 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");
+ fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.y*LENGTH_Angstrom/2,
region.y*LENGTH_Angstrom/2, "ylo yhi");
+ 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\n");
+ /* DO_SITE { */
+ /* fprintf(fp, "%7d%5d", n+1, site[n].type+1 ); // write atom-ID and
atom-type */
+ /* fprintf(fp, "%11.4f%11.4f%11.4f",
siteCoords[n].x*LENGTH_Angstrom, siteCoords[n].y*LENGTH_Angstrom,
siteCoords[n].z*LENGTH_Angstrom); */
+ /* fprintf(fp, "%7d", site[n].inLipid ); // write molecule-ID (this
is 0 for water sites) */
+ /* fprintf(fp, "%7.2f", charge[site[n].type]); // write charge (this
is 0.00 for water sites) */
+ /* fprintf(fp, "% 10.5f% 10.5f% 10.5f", // write dipole component
('dipole' style) */
+ /* 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, "%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"); */
+ /* } */
+ atomCount = 0;
+ DO_SITE {
+ atomCount++;
+ fprintf(fp, "%7d%5d", atomCount, site[n].type+1 ); // write atom-ID
and atom-type
+ fprintf(fp, "%11.4f%11.4f%11.4f", siteCoords[n].x*LENGTH_Angstrom,
siteCoords[n].y*LENGTH_Angstrom, siteCoords[n].z*LENGTH_Angstrom);
+ fprintf(fp, "%7d", site[n].inLipid ); // write molecule-ID (this is 0
for water sites)
+ fprintf(fp, "%7.2f", charge[site[n].type]); // write charge (this is
0.00 for water sites):
+ // write dipole component ('dipole' style):
+ if (site[n].type==GLYCEROL_TYPE || site[n].type==ESTER_TYPE) { //
generate two extra charged atoms to represent the dipole
+ fprintf(fp, "% 10.5f% 10.5f% 10.5f", 0.,0.,0.); // set dipole to
zero as it is going to be represented by two charges (below)
+ }
+ else { fprintf(fp, "% 10.5f% 10.5f% 10.5f", // write dipole components
+ 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, "%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");
+ if (site[n].type==GLYCEROL_TYPE || site[n].type==ESTER_TYPE) { //
generate two extra charged atoms to represent the dipole
+ disp = sigma[site[n].type] / 4; // displacement from dipole center
+ VSAdd( tailCoords, siteCoords[n], -disp, siteOrientVec[n] );
+ VSAdd( tipCoords, siteCoords[n], disp, siteOrientVec[n] );
+ atomCharge = dipole[site[n].type] * DIPOLE_D * D_IN_eA /
(2*disp*LENGTH_Angstrom);
+ // add "tail" charge atom:
+ atomCount++;
+ fprintf(fp, "%7d%5d", atomCount, site[n].type+1 ); // write atom-ID
and atom-type
+ fprintf(fp, "%11.4f%11.4f%11.4f", tailCoords.x*LENGTH_Angstrom,
tailCoords.y*LENGTH_Angstrom, tailCoords.z*LENGTH_Angstrom);
+ fprintf(fp, "%7d", site[n].inLipid ); // write molecule-ID
+ fprintf(fp, "%7.2f", -atomCharge); // write charge
+ fprintf(fp, "% 10.5f% 10.5f% 10.5f", 0.,0.,0.); // set dipole to zero
+ 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");
+ // add "tip" charge atom:
+ atomCount++;
+ fprintf(fp, "%7d%5d", atomCount, site[n].type+1 ); // write atom-ID
and atom-type
+ fprintf(fp, "%11.4f%11.4f%11.4f", tipCoords.x*LENGTH_Angstrom,
tipCoords.y*LENGTH_Angstrom, tipCoords.z*LENGTH_Angstrom);
+ fprintf(fp, "%7d", site[n].inLipid ); // write molecule-ID
+ fprintf(fp, "%7.2f", atomCharge); // write charge
+ fprintf(fp, "% 10.5f% 10.5f% 10.5f", 0.,0.,0.); // set dipole to zero
+ 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");
+ }
+ }
+ if (atomCount!=nAtoms) {printf("Error: atomCount!=nAtoms.\n"); exit(0);}
+ fprintf(fp, "\n");
+
+ fprintf(fp, "Bonds\n\n");
+ bondCount = 0;
+ DO_SITE {
+ id = n + (site[n].inLipid - 1) * 6; // conversion to LAMMPS atom-ID
+ 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 = id + 1;
+ atom2 = atom1 + 1;
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 1, atom1, atom2 );
+ } else { /* amine***phosphate (DOPE) */
+ atom1 = id + 1;
+ atom2 = atom1 + 1;
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 6, atom1, atom2 );
+ }
+ break;
+ case PHOSPHATE: // phosphate***glycerol
+ atom1 = id + 1;
+ atom2 = atom1 + 1;
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 2, atom1, atom2 );
+ break;
+ case GLYCEROL:
+ // GLYCEROL***ESTER_A:
+ atom1 = id + 1; // gly
+ // atom2 = atom1 + 1;
+ atom2 = atom1 + 3; // because of additional 2 sites used for
glycerol atom charges
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 3, atom1, atom2 );
+ // GLYCEROL***ESTER_B:
+ // atom2 = atom1 + 7;
+ atom2 = atom1 + 11; // because of additional 4 sites used for
glycerol and ester atom charges
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 3, atom1, atom2 );
+ // Gly COM *** dipole tail ('-' end)
+ atom2 = atom1 + 1;
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 9, atom1, atom2 );
+ // Gly COM *** dipole tip ('+' end)
+ atom2 = atom1 + 2;
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 9, atom1, atom2 );
+ // Gly dipole tip ('+' end) *** phosphate
+ atom1 = id; // phosphate (assumes phosphate'ID always precedes
glycerol's ID
+ atom2 = id + 3;
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 7, atom1, atom2 );
+ break;
+ case ESTER_A: // ESTER_A***TAIL_A1
+ case ESTER_B: // ESTER_B***TAIL_B1
+ if (site[n].lipidUnit == ESTER_A) estCom = id + 3;
+ else estCom = id + 5;
+ chain = estCom + 3;
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 4, estCom, chain );
+ estDipTail = estCom + 1;
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 8, estDipTail, chain );
+ bondCount++; // Est COM *** dipole tail:
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 9, estCom, estDipTail );
+ estDipTip = estCom + 2; // est tip
+ bondCount++; // Est COM *** dipole tip:
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 9, estCom, estDipTip );
+ 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: // TAIL_A4***TAIL_A5
+ case TAIL_B4: // TAIL_B4***TAIL_B5
+ atom1 = id + 7;
+ atom2 = atom1 + 1;
+ bondCount++;
+ fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 5, atom1, atom2 );
+ break;
+ default:
+ printf("Error: default case reached in WriteLammpsData.
Exiting...\n");
+ exit(0);
+ }
+ }
+ if (bondCount!=nBonds) {printf("Error:
bondCount=%d!=%d=nBonds.\n",bondCount,nBonds); exit(0);}
+ fprintf(fp, "\n");
+
+ fprintf(fp, "Angles\n\n");
+ angleCount = 0;
+ DO_SITE {
+ id = n + (site[n].inLipid - 1) * 6; // conversion to LAMMPS atom-ID
+ 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) */
+ angleCount++;
+ atom1 = id; // choline or amine
+ atom2 = id+1; // phosphate
+ atom3 = id+2; // glycerol
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 1, atom1, atom2, atom3
);
+ break;
+ case GLYCEROL:
+ angleCount++; // phos-gly-estA
+ atom1 = id; // phosphate
+ atom2 = id+1; // glycerol
+ atom3 = id+4; // ester_a
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 2, atom1, atom2, atom3
);
+ angleCount++; // phos-gly-estB
+ atom3 = id+12; // ester_b
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 2, atom1, atom2, atom3
);
+ angleCount++; /* dipTail-glyCOM-dipTip: */
+ atom1 = id+2; // dipTail
+ atom2 = id+1; // glycerol
+ atom3 = id+3; // dipTip
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 5, atom1, atom2, atom3
);
+ break;
+ case ESTER_A:
+ angleCount++; /* glycerol(3)-ester_A(4)-hydrocarbon_a1(5): */
+ atom1 = id; // gly
+ atom2 = id+3; // ester_a
+ atom3 = id+6; // tail_a1
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 3, atom1, atom2, atom3
);
+ angleCount++; /* dipTail-estCOM-dipTip: */
+ atom1 = id+4; // dipTail
+ atom2 = id+3; // ester_a
+ atom3 = id+5; // dipTip
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 5, atom1, atom2, atom3
);
+ break;
+ case TAIL_A1: /* ester_A(4)-hydrocarbon_a1(5)-hydrocarbon_a2(6) */
+ angleCount++;
+ atom1 = id+2; // estA
+ atom2 = id+5; // chain_a1
+ atom3 = id+6; // chain_a2
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 3, atom1, atom2, atom3
);
+ break;
+ case TAIL_B1: /* ester_B(4)-hydrocarbon_b1(5)-hydrocarbon_b2(6) */
+ angleCount++;
+ atom1 = id+4; // estB
+ atom2 = id+7; // chain_b1
+ atom3 = id+8; // chain_b2
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 3, atom1, atom2, atom3
);
+ break;
+ case TAIL_A3: /*
hydrocarbon_a2(6)-hydrocarbon_a3(7)-hydrocarbon_a4(8) */
+ case TAIL_A4: /*
hydrocarbon_a3(7)-hydrocarbon_a4(8)-hydrocarbon_a5(9) */
+ angleCount++;
+ atom1 = id+4;
+ atom2 = atom1+1;
+ atom3 = atom2+1;
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 3, atom1, atom2, atom3
);
+ break;
+ case TAIL_B3: /*
hydrocarbon_b2(6)-hydrocarbon_b3(7)-hydrocarbon_b4(8) */
+ case TAIL_B4: /*
hydrocarbon_b3(7)-hydrocarbon_b4(8)-hydrocarbon_b5(9) */
+ angleCount++;
+ atom1 = id+6;
+ atom2 = atom1+1;
+ atom3 = atom2+1;
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 3, atom1, atom2, atom3
);
+ break;
+ case ESTER_B: /* glycerol(3)-ester_B(10)-hydrocarbon_b1(11) -- SPECIAL
CASE */
+ angleCount++;
+ atom1 = id-6;
+ atom2 = id+5;
+ atom3 = id+8;
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 3, atom1, atom2, atom3
);
+ angleCount++; /* dipTail-estCOM-dipTip: */
+ atom1 = id+6; // dipTail
+ atom2 = id+5; // ester_b
+ atom3 = id+7; // dipTip
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 5, atom1, atom2, atom3
);
+ break;
+ case TAIL_A2: /*
hydrocarbon_a1(5)-hydrocarbon_a2(6)-hydrocarbon_a3(7) */
+ angleCount++;
+ atom1 = id+4;
+ atom2 = atom1+1;
+ atom3 = atom2+1;
+ // set cisUnsat angle type for DOPC and DOPE - case for DSPC to be
coded!
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 4, atom1, atom2, atom3
);
+ break;
+ case TAIL_B2: /*
hydrocarbon_b1(5)-hydrocarbon_b2(6)-hydrocarbon_b3(7) */
+ angleCount++;
+ atom1 = id+6;
+ atom2 = atom1+1;
+ atom3 = atom2+1;
+ // set cisUnsat angle type for DOPC and DOPE - case for DSPC to be
coded!
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 4, atom1, atom2, atom3
);
+ break;
+ default:
+ printf("Error in function WriteLammpsData - default case reached.
Exiting...\n"); exit(0);
+ }
+ }
+ if (angleCount!=nAngles) {printf("Error:
angleCount=%d!=%d=nAngles.\n",angleCount,nAngles); 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");
+}
=======================================
--- /trunk/lammpsScripts/lammpstrj2pdb.py Wed Jun 29 10:11:43 2011
+++ /trunk/lammpsScripts/lammpstrj2pdb.py Wed Jul 6 10:04:36 2011
@@ -17,23 +17,13 @@
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
+print "Setting dipole scaling factor assuming real units..."
+s=2; # scaling factor
#s=s/30; # rescale for reduced units - comment this line for real units
nLine=0; # line counter

@@ -46,41 +36,42 @@
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])
+ if len(words) == 9:
+ n=int(words[0]) # atom identifier
+ t=int(words[1]) # type identifier
+ m=float(words[2]) # molecule identifier
+ x=float(words[3])
+ y=float(words[4])
+ z=float(words[5])
+ mux=float(words[6])
+ muy=float(words[7])
+ muz=float(words[8])
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 ))
+ %(n,"OW WAT",m,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 ))
+ %(nPlus,"HW WAT",m,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 ))
+ %(n,"NC3 LIP",m,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 ))
+ %(n,"PO4 LIP",m,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 ))
+ %(n,"O LIP",m,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 ))
+ %(nPlus,"C LIP",m,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 ))
+ %(n,"CH2 LIP",m,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 ))
+ %(n,"NH3 LIP",m,x,y,z ))

print "Done scanning - end of script."
=======================================
--- /trunk/src/definitions.h Wed Jun 29 10:11:43 2011
+++ /trunk/src/definitions.h Wed Jul 6 10:04:36 2011
@@ -56,7 +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 D_IN_eA (0.2082) // 1 D = 0.2082 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/makefiles/Makefile_gcc-O3 Wed Jun 29 10:11:43 2011
+++ /trunk/src/makefiles/Makefile_gcc-O3 Wed Jul 6 10:04:36 2011
@@ -33,6 +33,7 @@

HEADERS = definitions.h dataStructs.h in_namelist.h
SOURCES = edp.c \
+ lammps.c \
latDiff.c \
lpp.c \
wpp.c \
@@ -59,6 +60,7 @@
topology.c \
utils.c
OBJECTS = edp.o \
+ lammps.o \
latDiff.o \
lpp.o \
wpp.o \
@@ -154,6 +156,8 @@

in_namelist.o: in_namelist.c dataStructs.h definitions.h

+lammps.o: lammps.c dataStructs.h definitions.h
+
latDiff.o: latDiff.c dataStructs.h definitions.h

lpp.o: lpp.c dataStructs.h definitions.h
=======================================
--- /trunk/src/makefiles/Makefile_icc-fast Wed Jun 29 10:11:43 2011
+++ /trunk/src/makefiles/Makefile_icc-fast Wed Jul 6 10:04:36 2011
@@ -33,6 +33,7 @@

in_namelist.o: in_namelist.c dataStructs.h definitions.h

+lammps.o: lammps.c dataStructs.h definitions.h
+
latDiff.o: latDiff.c dataStructs.h definitions.h

lpp.o: lpp.c dataStructs.h definitions.h
@@ -199,4 +203,3 @@
install:

uninstall:
-
=======================================
--- /trunk/src/output.c Wed Jun 29 10:11:43 2011
+++ /trunk/src/output.c Wed Jul 6 10:04:36 2011
@@ -19,20 +19,13 @@
extern const Atom *atom;
extern const LipidStruct dopcStruct;
extern const VecRTwo *xyLipidComTrue, xySoluteComTrue;
-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 real timeNow;
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, nDOPCsDSPCs, nDOPEs, nSolutes,
stepCount, nAtoms, zConstraint, nTypes;
+extern const int nSites, nLipids, nSolutes, stepCount;
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 );
@@ -41,7 +34,7 @@
void PrintBrahmsVersion()
{
printf("\n**********************\n");
- printf( "* BRAHMS - 28Jun2011 *\n");
+ printf( "* BRAHMS - 06Jul2011 *\n");
printf( "**********************\n\n");
fflush(stdout);
}
@@ -376,238 +369,3 @@
{
fprintf( fPtrAV, "%.3f %.2f\n", xyArea.val, boxVolume.val ); /* [nm^2],
[nm^3] */
}
-
-/*******************************************************************************************************
- * WriteLammpsData -- writes "data." file to be read into the LAMMPS
program through read_data command *
-
*******************************************************************************************************/
-void WriteLammpsData( FILE *fp )
-{
- 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");
- fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.y*LENGTH_Angstrom/2,
region.y*LENGTH_Angstrom/2, "ylo yhi");
- 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\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, "%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, "%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, "% 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, "%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");
-}

Reply all

Reply to author

Forward

0 new messages