[brahms-md] r53 committed - Working on better representation of lipid dipoles for LAMMPS.
brah...@googlecode.com
Author: orsimario
Date: Thu Aug 18 10:23:37 2011
Log: Working on better representation of lipid dipoles for LAMMPS.
http://code.google.com/p/brahms-md/source/detail?r=53
Modified:
/trunk/src/forceBdMain.c
/trunk/src/forceNbdMain.c
/trunk/src/lammps.c
/trunk/src/output.c
=======================================
--- /trunk/src/forceBdMain.c Wed Jun 29 10:11:43 2011
+++ /trunk/src/forceBdMain.c Thu Aug 18 10:23:37 2011
@@ -1,4 +1,4 @@
-/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
+/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 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.
=======================================
--- /trunk/src/forceNbdMain.c Wed Jul 20 07:59:24 2011
+++ /trunk/src/forceNbdMain.c Thu Aug 18 10:23:37 2011
@@ -62,8 +62,8 @@
if ( interactionType[ n ] == WAT_WAT ) { /* water-water interaction
(special case) */
if ( rr < rrCutWatWat ) {
//
CalcDipDipInteraction(&ijPotEn,&ijForceTot,&ijTorqueTot,&jiTorqueTot,watDipoleSq,siteOrientVec[i],siteOrientVec[j],rij,rr);
-// CalcDipDipInteraction_CS( &ijPotEn, &ijForceTot, &ijTorqueTot,
&jiTorqueTot, watDipoleSq, siteOrientVec[i], siteOrientVec[j], rCutWatWat,
rij, rr );
- CalcDipDipInteraction_SF( &ijPotEn, &ijForceTot, &ijTorqueTot,
&jiTorqueTot, watDipoleSq, siteOrientVec[i], siteOrientVec[j], rCutWatWat,
rij, rr );
+ CalcDipDipInteraction_CS( &ijPotEn, &ijForceTot, &ijTorqueTot,
&jiTorqueTot, watDipoleSq, siteOrientVec[i], siteOrientVec[j], rCutWatWat,
rij, rr );
+ // CalcDipDipInteraction_SF( &ijPotEn, &ijForceTot, &ijTorqueTot,
&jiTorqueTot, watDipoleSq, siteOrientVec[i], siteOrientVec[j], rCutWatWat,
rij, rr );
potEnWatSum += ijPotEn;
// CalcLennJonesInteraction( &ijPotEn, &ijForce, sigWatSq, epsWat,
rr, rij );
CalcLennJonesInteraction_SF( &ijPotEn, &ijForce, sigWatSq, epsWat, rr,
rrCutWatWat, rij );
=======================================
--- /trunk/src/lammps.c Wed Jul 20 07:59:24 2011
+++ /trunk/src/lammps.c Thu Aug 18 10:23:37 2011
@@ -15,32 +15,29 @@
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 real springRigidity, angleRigidity, orientRigidity,
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; // types defined below:
+static const int nAtomsPerLipid = 15; // assuming 15-site lipid models
+static const int nHarmonicBondsPerLipid = 14; // assuming 15-site lipid
models
+static const int nDipoleBondsPerLipid = 3; // assuming 3 dipoles to
restrain
+
+static const int nBondTypes = 6; // 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 DIPOLE = 6;
+static const int AMINE_PHOS = 7;
+
+static const int nAnglesPerLipid = 13; // assuming 15-site lipid models: 13
+static const int nAngleTypes = 4; // 1_chol(amine)-phos-gly,
2_phos-gly-est, 3_PI_sat, 4_cisUnsat
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;
/*******************************************************************************************************
@@ -49,13 +46,11 @@
void WriteLammpsData( FILE *fp )
{
int n, id, atomCount, bondCount, angleCount, atom1, atom2, atom3;
- int idPhos, idGly, idEstA, idEstB, idGlyDipTip, idGlyDipTail, idEst,
estDipTail, estDipTip, idChain1, idChainI, idChainJ;
+ int idGly, idEstA, idEstB, idEst, idChain1, idChainI, idChainJ;
int nAtoms = nWaters + nAtomsPerLipid*nLipids; // #atoms for LAMMPS
topology
- int nBonds = nBondsPerLipid*nLipids;
+ int nBonds = (nHarmonicBondsPerLipid+nDipoleBondsPerLipid)*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
@@ -78,61 +73,22 @@
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
+ if ( site[n].type==AMINE_TYPE )
+ fprintf(fp, "%7d%5d", atomCount, 2 ); // write atom-ID and atom-type
for amine - TO SORT!
+ else 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, 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
(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 );
- }
+ fprintf(fp, "% 10.5f% 10.5f% 10.5f", // write dipole components
(non-zero only for water, glycerol and ester 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 );
// 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( 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 "tip" charge atom:
- atomCount++;
- 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 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 "tail" charge atom:
- atomCount++;
- 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, "% 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");
@@ -140,7 +96,7 @@
fprintf(fp, "Bonds\n\n");
bondCount = 0;
DO_SITE {
- id = n + (site[n].inLipid - 1) * 6; // conversion to LAMMPS atom-ID
+ id = n + 1 + (site[n].inLipid - 1); // 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
@@ -148,50 +104,44 @@
break;
case CHOLINE_OR_AMINE: /* choline***phosphate (DOPC) or
amine***phosphate (DOPE) */
if (site[n].type == CHOLINE_TYPE) { /* choline***phosphate (DOPC) */
- atom1 = id + 1;
+ atom1 = id;
atom2 = atom1 + 1;
fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHOL_PHOS, atom1, atom2 );
} else { /* amine***phosphate (DOPE) */
- atom1 = id + 1;
+ atom1 = id;
atom2 = atom1 + 1;
- // fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, AMINE_PHOS, atom1, atom2 );
TO BE SORTED!
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, AMINE_PHOS, atom1, atom2 );
}
break;
case PHOSPHATE: // phosphate***glycerol
- atom1 = id + 1;
- atom2 = atom1 + 1;
+ atom1 = id; // phos
+ atom2 = atom1 + 1; // gly
fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, PHOS_GLY, atom1, atom2 );
+ // "dipole" bond:
+ // fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, DIPOLE, atom1,
atom2 ); // gly dipole -- phos
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, DIPOLE, atom1-1, atom2 );
// gly dipole -- chol
break;
case GLYCEROL:
- idGly = id + 1; // gly
- idEstA = idGly + 3; // because of additional 2 sites used for
glycerol atom charges
+ idGly = id; // gly
+ idEstA = idGly + 1;
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
+ idEstB = idGly + 7; // because there are 7 sites between glycerol
and esterB in the topological definition of a lipid molecule
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) idEst = id + 3;
- else idEst = id + 5;
- idChain1 = idEst + 3;
+ idEst = id;
+ idChain1 = idEst + 1;
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
+ // "dipole" bond:
+ // fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, DIPOLE, idChain1,
idEst ); // chain1 -- estDip
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, DIPOLE, idChain1+1, idEst
); // chain2 -- estDip
break;
case TAIL_A1: // TAIL_A1***TAIL_A2
case TAIL_A2: // TAIL_A2***TAIL_A3
case TAIL_A3: // TAIL_A3***TAIL_A4
case TAIL_A4: // TAIL_A4***TAIL_A5
- idChainI = id + 5;
+ idChainI = id;
idChainJ = idChainI + 1;
fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHAIN_CHAIN, idChainI,
idChainJ );
break;
@@ -199,7 +149,7 @@
case TAIL_B2: // TAIL_B2***TAIL_B3
case TAIL_B3: // TAIL_B3***TAIL_B4
case TAIL_B4: // TAIL_B4***TAIL_B5
- idChainI = id + 7;
+ idChainI = id;
idChainJ = idChainI + 1;
fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHAIN_CHAIN, idChainI,
idChainJ );
break;
@@ -214,7 +164,7 @@
fprintf(fp, "Angles\n\n");
angleCount = 0;
DO_SITE {
- id = n + (site[n].inLipid - 1) * 6; // conversion to LAMMPS atom-ID
+ id = n + (site[n].inLipid - 1); // 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:
@@ -230,72 +180,34 @@
case GLYCEROL:
atom1 = id; // phosphate
atom2 = id+1; // glycerol
- atom3 = id+4; // ester_a
+ atom3 = id+2; // ester_a
fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2,
atom3 ); // phos-gly-estA
- atom3 = id+12; // ester_b
+ atom3 = id+8; // ester_b
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, _piDipCharges_,
atom1, atom2, atom3 ); /* dipTail-glyCOM-dipTip: */
break;
- case ESTER_A:
- atom1 = id; // gly
- atom2 = id+3; // ester_a
- atom3 = id+6; // tail_a1
- 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, _piDipCharges_,
atom1, atom2, atom3 ); /* dipTail-estCOM-dipTip: */
- break;
- 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, _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) */
- atom1 = id+4; // estB
- atom2 = id+7; // chain_b1
- atom3 = id+8; // chain_b2
- fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, 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_B1: /* ester_B(4)-hydrocarbon_b1(5)-hydrocarbon_b2(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) */
- atom1 = id+4;
- atom2 = atom1+1;
- atom3 = atom2+1;
- 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) */
- atom1 = id+6;
- atom2 = atom1+1;
- atom3 = atom2+1;
+ atom1 = id;
+ atom2 = id+1;
+ atom3 = id+2;
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 */
atom1 = id-6;
- atom2 = id+5;
- atom3 = id+8;
+ atom2 = id+1;
+ atom3 = id+2;
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, _piDipCharges_,
atom1, atom2, atom3 ); /* dipTail-estCOM-dipTip: */
break;
case TAIL_A2: /*
hydrocarbon_a1(5)-hydrocarbon_a2(6)-hydrocarbon_a3(7) */
- 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, _cisUnsat_, atom1,
atom2, atom3 );
- break;
case TAIL_B2: /*
hydrocarbon_b1(5)-hydrocarbon_b2(6)-hydrocarbon_b3(7) */
- atom1 = id+6;
- atom2 = atom1+1;
- atom3 = atom2+1;
+ atom1 = id;
+ atom2 = id+1;
+ atom3 = id+2;
// set cisUnsat angle type for DOPC and DOPE - case for DSPC to be
coded!
fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _cisUnsat_, atom1,
atom2, atom3 );
break;
@@ -307,30 +219,31 @@
}
/**********************************************************************************************
- * WriteLammpsForcefield -- writes "forcefield." file to be included into
LAMMPS input script *
+ * WriteLammpsForcefield -- writes "forcefield" file to include 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 orientRigidity_LAMMPS_units = orientRigidity * 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, "pair_style dipole/sf 12\n");
+ fprintf(fp, "bond_style hybrid harmonic dipole\n");
fprintf(fp, "angle_style cosine/squared\n");
fprintf(fp, "\n");
- for ( n = 0; n < nTypes; n++ ) { // write masses
+ for ( n = 0; n < nTypes-1; 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++ ) {
+ for ( i = 0; i < nTypes-1; i++ ) {
+ for ( j = 0; j < nTypes-1; 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);
@@ -343,16 +256,13 @@
// 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", 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 harmonic%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 harmonic%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 harmonic%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 harmonic%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 harmonic%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 dipole %7.3f%7.3f\n", DIPOLE,
orientRigidity_LAMMPS_units/2, 0.); // '0' is irrelevant at the moment, as
gamma_reference = 0
// 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); */
@@ -367,7 +277,6 @@
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/output.c Wed Jul 20 07:59:24 2011
+++ /trunk/src/output.c Thu Aug 18 10:23:37 2011
@@ -34,7 +34,7 @@
void PrintBrahmsVersion()
{
printf("\n**********************\n");
- printf( "* BRAHMS - 20Jul2011 *\n");
+ printf( "* BRAHMS - 18Aug2011 *\n");
printf( "**********************\n\n");
fflush(stdout);
}