[brahms-md] r51 committed - Fixed bug in the preparation of bond topology for LAMMPS, made various...
brah...@googlecode.com
Author: orsimario
Date: Fri Jul 8 09:15:59 2011
Log: Fixed bug in the preparation of bond topology for LAMMPS, made
various cosmetic changes for improved clarity
http://code.google.com/p/brahms-md/source/detail?r=51
Modified:
/trunk/lammpsScripts/lammpstrj2pdb.py
/trunk/src/dataStructs.h
/trunk/src/lammps.c
/trunk/src/topology.c
=======================================
--- /trunk/lammpsScripts/lammpstrj2pdb.py Wed Jul 6 10:04:36 2011
+++ /trunk/lammpsScripts/lammpstrj2pdb.py Fri Jul 8 09:15:59 2011
@@ -64,14 +64,14 @@
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 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 LIP",m,xh,yh,zh ))
+# 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 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 LIP",m,x,y,z ))
- elif t == 7: # amine type
+ elif t == 7: # charge type
outFile.write('ATOM%7d%9s%6d%12.3f%8.3f%8.3f\n'
- %(n,"NH3 LIP",m,x,y,z ))
+ %(n,"C LIP",m,x,y,z ))
print "Done scanning - end of script."
=======================================
--- /trunk/src/dataStructs.h Tue Jun 28 10:57:32 2011
+++ /trunk/src/dataStructs.h Fri Jul 8 09:15:59 2011
@@ -39,7 +39,8 @@
ESTER_TYPE = 4,
TAIL_TYPE = 5,
AMINE_TYPE = 6,
- SOLUTE_TYPE = 7,
+ CHARGE_TYPE = 7, // for the additional charges (with zero LJ
parameters) used to represent est and gly dipoles in LAMMPS
+ SOLUTE_TYPE = 8
} SiteType;
typedef enum { WAT_WAT = 0,
=======================================
--- /trunk/src/lammps.c Wed Jul 6 10:04:36 2011
+++ /trunk/src/lammps.c Fri Jul 8 09:15:59 2011
@@ -21,9 +21,27 @@
// 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 nBondTypes = 10; // types defined below:
+static const int CHOL_PHOS = 1;
+static const int PHOS_GLY = 2;
+static const int GLY_EST =3;
+static const int EST_CHAIN = 4;
+static const int CHAIN_CHAIN = 5;
+//static const int AMINE_PHOS = 6;
+static const int GLYPLUS_PHOS = 6;
+static const int ESTPLUS_CHAIN = 7;
+static const int GLYCHARGE_COM = 8;
+static const int ESTCHARGE_COM = 9;
+
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
+static const int _phos_ = 1;
+static const int _gly_ = 2;
+static const int _piSat_ = 3;
+static const int _cisUnsat_ = 4;
+static const int _piDipCharges_ = 5;
+
/*******************************************************************************************************
* WriteLammpsData -- writes "data." file to be read into the LAMMPS
program through read_data command *
@@ -31,13 +49,17 @@
void WriteLammpsData( FILE *fp )
{
int n, id, atomCount, bondCount, angleCount, atom1, atom2, atom3;
- int estCom, estDipTail, estDipTip, chain;
+ int idPhos, idGly, idEstA, idEstB, idGlyDipTip, idGlyDipTail, idEst,
estDipTail, estDipTip, idChain1, idChainI, idChainJ;
int nAtoms = nWaters + nAtomsPerLipid*nLipids; // #atoms for LAMMPS
topology
int nBonds = nBondsPerLipid*nLipids;
int nAngles = nAnglesPerLipid*nLipids;
+ int nAtomTypes;
real disp, atomCharge;
VecR tailCoords, tipCoords;
+ if ( nLipids == nDOPCsDSPCs || nLipids == nDOPEs ) nAtomTypes =
nTypes-1; // bilayer containing only one lipid species
+ else nAtomTypes = nTypes; // mixed bilayer
+
fprintf(fp, "LAMMPS 'data.' description - created with Brahms
(www.personal.soton.ac.uk/orsi/brahms)\n");
fprintf(fp, "\n");
@@ -45,7 +67,7 @@
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 atom types\n", nAtomTypes);
fprintf(fp, "%9d bond types\n", nBondTypes);
fprintf(fp, "%9d angle types\n", nAngleTypes);
fprintf(fp, "\n");
@@ -74,13 +96,12 @@
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):
+ fprintf(fp, "%7d", site[n].inLipid ); // write molecule-ID (this is 0
for water sites, 1 for lipid-1, 2 for lipid-2, etc.)
+ fprintf(fp, "%7.2f", charge[site[n].type]); // write charge (this is
0.00 for point-dipole sites):
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
+ else { fprintf(fp, "% 10.5f% 10.5f% 10.5f", // write dipole components
(non-zero only for water sites)
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 );
@@ -90,24 +111,24 @@
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] );
+ VSAdd( tailCoords, siteCoords[n], -disp, siteOrientVec[n] );
atomCharge = dipole[site[n].type] * DIPOLE_D * D_IN_eA /
(2*disp*LENGTH_Angstrom);
- // add "tail" charge atom:
+ // 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", tailCoords.x*LENGTH_Angstrom,
tailCoords.y*LENGTH_Angstrom, tailCoords.z*LENGTH_Angstrom);
+ fprintf(fp, "%7d%5d", atomCount, CHARGE_TYPE ); // 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, "%7.2f", atomCharge); // write charge
+ fprintf(fp, "% 10.5f% 10.5f% 10.5f", 0.,0.,0.); // set point-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:
+ // 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", tipCoords.x*LENGTH_Angstrom,
tipCoords.y*LENGTH_Angstrom, tipCoords.z*LENGTH_Angstrom);
+ fprintf(fp, "%7d%5d", atomCount, CHARGE_TYPE ); // 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, "%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");
@@ -129,75 +150,58 @@
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 );
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHOL_PHOS, atom1, atom2 );
} else { /* amine***phosphate (DOPE) */
atom1 = id + 1;
atom2 = atom1 + 1;
- bondCount++;
- fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 6, atom1, atom2 );
+ // fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, AMINE_PHOS, atom1, atom2 );
TO BE SORTED!
}
break;
case PHOSPHATE: // phosphate***glycerol
atom1 = id + 1;
atom2 = atom1 + 1;
- bondCount++;
- fprintf(fp, "%7d%5d%8d%8d\n", bondCount, 2, atom1, atom2 );
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, PHOS_GLY, 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 );
+ idGly = id + 1; // gly
+ idEstA = idGly + 3; // because of additional 2 sites used for
glycerol atom charges
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLY_EST, idGly, idEstA );
// GLYCEROL***ESTER_A:
+ idEstB = idGly + 11; // because of additional 4 sites used for
glycerol and ester atom charges
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLY_EST, idGly, idEstB );
// GLYCEROL***ESTER_B:
+ idGlyDipTip = idGly + 1;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLYCHARGE_COM, idGly,
idGlyDipTip ); // Gly COM *** dipole tip ('+' end)
+ idGlyDipTail = idGly + 2;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLYCHARGE_COM, idGly,
idGlyDipTail ); // Gly COM *** dipole tail ('-' end)
+ idPhos = id; // phosphate (assumes phosphate'ID always precedes
glycerol's ID
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLYPLUS_PHOS, idPhos,
idGlyDipTip ); // Gly dipole tip ('+' end) *** phosphate
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 );
+ if (site[n].lipidUnit == ESTER_A) idEst = id + 3;
+ else idEst = id + 5;
+ idChain1 = idEst + 3;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, EST_CHAIN, idEst,
idChain1 );
+ estDipTail = idEst + 2;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, ESTPLUS_CHAIN,
estDipTail, idChain1 );
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, ESTCHARGE_COM, idEst,
estDipTail ); // Est COM *** dipole tail
+ estDipTip = idEst + 1; // est tip
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, ESTCHARGE_COM, idEst,
estDipTip ); // Est COM *** dipole tip
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
+ idChainI = id + 5;
+ idChainJ = idChainI + 1;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHAIN_CHAIN, idChainI,
idChainJ );
+ break;
+ case TAIL_B1: // TAIL_B1***TAIL_B2
+ case TAIL_B2: // TAIL_B2***TAIL_B3
+ case TAIL_B3: // TAIL_B3***TAIL_B4
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 );
+ idChainI = id + 7;
+ idChainJ = idChainI + 1;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHAIN_CHAIN, idChainI,
idChainJ );
break;
default:
printf("Error: default case reached in WriteLammpsData.
Exiting...\n");
@@ -218,96 +222,82 @@
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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _phos_, 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
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2,
atom3 ); // phos-gly-estA
atom3 = id+12; // ester_b
- fprintf(fp, "%7d%5d%8d%8d%8d\n", angleCount, 2, atom1, atom2, atom3
);
- angleCount++; /* dipTail-glyCOM-dipTip: */
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2,
atom3 ); // phos-gly-estB
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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piDipCharges_,
atom1, atom2, atom3 ); /* dipTail-glyCOM-dipTip: */
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: */
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1,
atom2, atom3 ); /* glycerol(3)-ester_A(4)-hydrocarbon_a1(5): */
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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piDipCharges_,
atom1, atom2, atom3 ); /* dipTail-estCOM-dipTip: */
break;
- case TAIL_A1: /* ester_A(4)-hydrocarbon_a1(5)-hydrocarbon_a2(6) */
- angleCount++;
+ case TAIL_A1:
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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1,
atom2, atom3 ); /* ester_A(4)-hydrocarbon_a1(5)-hydrocarbon_a2(6) */
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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, 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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, 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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, 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: */
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1,
atom2, atom3 );
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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piDipCharges_,
atom1, atom2, atom3 ); /* dipTail-estCOM-dipTip: */
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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _cisUnsat_, 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
);
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _cisUnsat_, atom1,
atom2, atom3 );
break;
default:
printf("Error in function WriteLammpsData - default case reached.
Exiting...\n"); exit(0);
@@ -352,16 +342,16 @@
// 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", CHOL_PHOS,
springRigidity_LAMMPS_units/2,
0.9*sig[CHOLINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", PHOS_GLY,
springRigidity_LAMMPS_units/2,
0.9*sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLY_EST,
springRigidity_LAMMPS_units/2,
0.9*sig[GLYCEROL_TYPE][ESTER_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", EST_CHAIN,
springRigidity_LAMMPS_units/2,
0.9*sig[ESTER_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", CHAIN_CHAIN,
springRigidity_LAMMPS_units/2,
0.9*sig[TAIL_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
+ // fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", AMINE_PHOS,
springRigidity_LAMMPS_units/2,
0.9*sig[PHOSPHATE_TYPE][AMINE_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLYPLUS_PHOS,
springRigidity_LAMMPS_units/2,
(sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]-sigma[GLYCEROL_TYPE]/4)*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", ESTPLUS_CHAIN,
springRigidity_LAMMPS_units/2,
(sig[ESTER_TYPE][TAIL_TYPE]-sigma[ESTER_TYPE]/4)*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLYCHARGE_COM,
springRigidity_LAMMPS_units/2, sigma[GLYCEROL_TYPE]/4*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", ESTCHARGE_COM,
springRigidity_LAMMPS_units/2, sigma[ESTER_TYPE]/4*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); */
@@ -372,10 +362,10 @@
// 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, "angle_coeff%5d%8.4f%8.2f\n", _phos_,
angleRigidity_LAMMPS_units/2, refAngle_CholPhosGly*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _gly_,
angleRigidity_LAMMPS_units/2, refAngle_PhosGlyEst*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _piSat_,
angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _cisUnsat_,
angleRigidity_LAMMPS_units/2, refAngle_cisUnsat*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _piDipCharges_,
angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
fprintf(fp, "\n");
}
=======================================
--- /trunk/src/topology.c Tue Jun 28 10:57:32 2011
+++ /trunk/src/topology.c Fri Jul 8 09:15:59 2011
@@ -109,7 +109,12 @@
charge[ type ] = 0;
dipole[ type ] = 0;
break;
- case SOLUTE_TYPE: // temporary fix to avoid compilation error
+ case CHARGE_TYPE:
+ sigma[ type ] = 0;
+ epsilon[ type ] = 0;
+ electroMoment[ type ] = 0;
+ charge[ type ] = 0;
+ dipole[ type ] = 0;
break;
default: if ( !zConstraint ) printf( "\nERROR in LoadTopolgy() - got
to the switch default case!\n\n" );
}