[brahms-md] r59 committed - Removed obsolete variables 'insertSolute', 'sigCharge' and 'epsCharge'...
brah...@googlecode.com
Author: orsimario
Date: Wed Sep 28 10:00:59 2011
Log: Removed obsolete variables 'insertSolute', 'sigCharge'
and 'epsCharge'. Introduced input parameter 'zReplicas', to request system
replication along the z-dimension; this is useful to generate multilamellar
systems. Improvements to 'lammps.c': ability to output z-replicated
systems, handling of new LAMMPS 'angle_dipole' style.
http://code.google.com/p/brahms-md/source/detail?r=59
Added:
/trunk/src/depot/lammps_old_dipole2charges.c
Modified:
/trunk/src/epp.c
/trunk/src/in_namelist.c
/trunk/src/in_namelist.h
/trunk/src/lammps.c
/trunk/src/main.c
/trunk/src/output.c
/trunk/src/startUp.c
/trunk/src/stepMain.c
=======================================
--- /dev/null
+++ /trunk/src/depot/lammps_old_dipole2charges.c Wed Sep 28 10:00:59 2011
@@ -0,0 +1,375 @@
+/* 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; // 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 = 1;
+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 *
+
*******************************************************************************************************/
+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 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");
+
+ 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", nAtomTypes);
+ 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++;
+ if ( site[n].type==AMINE_TYPE )
+ fprintf(fp, "%7d%5d", atomCount, 2 ); // write atom-ID and atom-type
foe 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 );
+ }
+ // 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");
+
+ 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;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHOL_PHOS, atom1, atom2 );
+ } else { /* amine***phosphate (DOPE) */
+ atom1 = id + 1;
+ atom2 = atom1 + 1;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, AMINE_PHOS, atom1, atom2 );
+ }
+ break;
+ case PHOSPHATE: // phosphate***glycerol
+ atom1 = id + 1;
+ atom2 = atom1 + 1;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, PHOS_GLY, atom1, atom2 );
+ break;
+ case GLYCEROL:
+ 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) 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_A2: // TAIL_A2***TAIL_A3
+ case TAIL_A3: // TAIL_A3***TAIL_A4
+ 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
+ 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");
+ 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) */
+ atom1 = id; // choline or amine
+ atom2 = id+1; // phosphate
+ atom3 = id+2; // glycerol
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _phos_, atom1, atom2,
atom3 );
+ break;
+ case GLYCEROL:
+ atom1 = id; // phosphate
+ atom2 = id+1; // glycerol
+ atom3 = id+4; // ester_a
+ 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, _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 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;
+ 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;
+ 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;
+ // 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;
+ 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");
+
+ if (nSites!=nWaters) { // to prevent segmentation fault when simulating
pure water systems
+ // 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%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); */
+ /* 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", _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/epp.c Thu Jan 6 06:46:39 2011
+++ /trunk/src/epp.c Wed Sep 28 10:00:59 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/in_namelist.c Wed Jul 20 07:59:24 2011
+++ /trunk/src/in_namelist.c Wed Sep 28 10:00:59 2011
@@ -24,7 +24,7 @@
real springRigidity;
int diffusion, /* logical: 1-> evaluate diffusion, 0-> don't do it */
writeLipLatMotion, stepDiffuse, latDiff, epp, stepEpp,
sizeHistEpp,limitWatRdf,sizeHistWatRdf;
-int limitEdp, limitLpp, limitEpp, limitWpp, limitDiff, centerInputStruct,
randSeed, insertSolute;
+int limitEdp, limitLpp, limitEpp, limitWpp, limitDiff, centerInputStruct,
randSeed, zReplicas;
int applyThermostat, applyBarostat, flexBox, stepNebr;
int adjustRegion; // input parameter; logical, used to require
modification of the region edges loaded from previous run
int resetTime; // input parameter; logical, used to require resetting
of stepCount and timeNow variables
@@ -49,8 +49,8 @@
extern VecR region, regionAdjusted, solInertia;
extern VecI initUcell, initHalfCellWat;
extern real scaleRefSpringLength;
-extern real sigChol, sigAmine, sigPhos, sigEst, sigGly, sigWat, sigCharge;
-extern real epsChol, epsAmine, epsPhos, epsGly, epsWat, epsCharge;
+extern real sigChol, sigAmine, sigPhos, sigEst, sigGly, sigWat;
+extern real epsChol, epsAmine, epsPhos, epsGly, epsWat;
extern real hBondAmineGly, hBondWatWat;
extern real hBondWatEst, hBondWatPhos, hBondWatGly;
extern real hBondAmineWat, hBondAminePhos, hBondAmineEst;
@@ -133,7 +133,7 @@
NameI( centerInputStruct),
NameI( reCenterBilayer ),
NameI( zConstraint ),
- NameI( insertSolute ),
+ NameI( zReplicas ),
};
NameList lipidForcefieldList[] = {
@@ -144,14 +144,12 @@
NameR( sigGly ),
NameR( sigEst ),
NameR( sigTail ),
- NameR( sigCharge ),
NameR( epsChol ),
NameR( epsAmine ),
NameR( epsPhos ),
NameR( epsGly ),
NameR( epsEst ),
NameR( epsTail ),
- NameR( epsCharge ),
NameR( hBondWatPhos ),
NameR( hBondWatGly ),
NameR( hBondWatEst ),
=======================================
--- /trunk/src/in_namelist.h Thu Jan 6 06:46:39 2011
+++ /trunk/src/in_namelist.h Wed Sep 28 10:00:59 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/lammps.c Sat Aug 20 05:29:45 2011
+++ /trunk/src/lammps.c Wed Sep 28 10:00:59 2011
@@ -14,30 +14,32 @@
extern const VecR region, *siteCoords, *siteOrientVec;
extern const Site *site;
-extern const real timeNow, *mass, *charge, *dipole, *sigma, **sig, **eps;
+extern const real timeNow, *mass, *charge, *dipole, *sigma, **sig, **eps,
massAmine;
extern const real springRigidity, angleRigidity, orientRigidity,
refAngle_CholPhosGly, refAngle_PhosGlyEst, refAngle_cisUnsat;
extern const int nSites, nLipids, nDOPCsDSPCs, nDOPEs, nTypes, nWaters;
+extern int zReplicas;
// local variables:
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 nBondTypes = 5; // types defined below:
+static const int X_PHOS = 1;
+//static const int CHOL_PHOS = 1;
+//static const int AMINE_PHOS = 7;
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 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 nCosineAnglesPerLipid = 13; // assuming 15-site lipid
models: 13
+static const int nDipoleAnglesPerLipid = 3; // assuming 3 dipoles to
restrain
+static const int nAngleTypes = 5; // 1_chol(amine)-phos-gly,
2_phos-gly-est, 3_PI_sat, 4_cisUnsat, 5_dipoles
static const int _phos_ = 1;
static const int _gly_ = 2;
static const int _piSat_ = 3;
static const int _cisUnsat_ = 4;
+static const int DIPOLE = 5;
/*******************************************************************************************************
@@ -45,16 +47,22 @@
*******************************************************************************************************/
void WriteLammpsData( FILE *fp )
{
- int n, id, atomCount, bondCount, angleCount, atom1, atom2, atom3;
+ int n, id, atomCount, bondCount, angleCount, atom1, atom2, atom3, zRep;
int idGly, idEstA, idEstB, idEst, idChain1, idChainI, idChainJ;
int nAtoms = nWaters + nAtomsPerLipid*nLipids; // #atoms for LAMMPS
topology
- int nBonds = (nHarmonicBondsPerLipid+nDipoleBondsPerLipid)*nLipids;
- int nAngles = nAnglesPerLipid*nLipids;
+ int nBonds = nHarmonicBondsPerLipid*nLipids;
+ int nAngles = (nCosineAnglesPerLipid+nDipoleAnglesPerLipid)*nLipids;
int nAtomTypes;
if ( nLipids == nDOPCsDSPCs || nLipids == nDOPEs ) nAtomTypes =
nTypes-1; // bilayer containing only one lipid species
else nAtomTypes = nTypes; // mixed bilayer
+ if ( zReplicas==0 ) zReplicas = 1; // set zReplicas to default value in
case it wasn't set in 'brahms.md'
+
+ nAtoms *= zReplicas;
+ nBonds *= zReplicas;
+ nAngles *= zReplicas;
+
fprintf(fp, "LAMMPS 'data.' description - created with Brahms
(www.personal.soton.ac.uk/orsi/brahms)\n");
fprintf(fp, "\n");
@@ -69,146 +77,148 @@
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, "% 13.7G% 13.7G%13s\n", -region.z*LENGTH_Angstrom/2,
region.z*(zReplicas-0.5)*LENGTH_Angstrom, "zlo zhi");
fprintf(fp, "\n");
fprintf(fp, "Atoms\n\n");
atomCount = 0;
- DO_SITE {
- atomCount++;
- 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):
- 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 (atomCount!=nAtoms) {printf("Error: atomCount!=nAtoms.\n"); exit(0);}
+ for ( zRep = 0; zRep < zReplicas; zRep++ ) {
+ DO_SITE {
+ atomCount++;
+ if ( site[n].type==AMINE_TYPE )
+ fprintf(fp, "%7d%5d", atomCount, 2 );
+ 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+zRep*region.z)*LENGTH_Angstrom);
+ if ( site[n].inLipid == 0 ) // water
+ fprintf(fp, "%7d", site[n].inLipid); // write molecule-ID (this is 0 for
water sites, 1 for lipid-1, 2 for lipid-2, etc.)
+ else fprintf(fp, "%7d", site[n].inLipid+(zRep*nLipids) ); // 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):
+ 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 (atomCount!=nAtoms) {printf("Error: atomCount = %d != %d =
nAtoms.\n", atomCount, nAtoms); exit(0);}
fprintf(fp, "\n");
fprintf(fp, "Bonds\n\n");
bondCount = 0;
- DO_SITE {
- id = n + 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
- 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) */
- fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHOL_PHOS, id, id+1 );
- else /* amine***phosphate (DOPE) */
- fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, AMINE_PHOS, id, id+1 );
- break;
- case PHOSPHATE: // phosphate***glycerol
- 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; // gly
- idEstA = idGly + 1;
- fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLY_EST, idGly, idEstA );
// GLYCEROL***ESTER_A:
- 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:
- break;
- case ESTER_A: // ESTER_A***TAIL_A1
- case ESTER_B: // ESTER_B***TAIL_B1
- idEst = id;
- idChain1 = idEst + 1;
- fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, EST_CHAIN, idEst,
idChain1 );
- // "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;
- 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
- idChainI = id;
- 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");
- exit(0);
- }
+ for ( zRep = 0; zRep < zReplicas; zRep++ ) {
+ DO_SITE {
+ id = (n + 1) + zRep * ( nAtoms / zReplicas ); // 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) */
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, X_PHOS, id, id+1 );
+ break;
+ case PHOSPHATE: // phosphate***glycerol
+ atom1 = id; // phos
+ atom2 = atom1 + 1; // gly
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, PHOS_GLY, atom1, atom2 );
+ break;
+ case GLYCEROL:
+ idGly = id; // gly
+ idEstA = idGly + 1;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLY_EST, idGly, idEstA ); //
GLYCEROL***ESTER_A:
+ 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:
+ break;
+ case ESTER_A: // ESTER_A***TAIL_A1
+ case ESTER_B: // ESTER_B***TAIL_B1
+ idEst = id;
+ idChain1 = idEst + 1;
+ fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, EST_CHAIN, idEst, idChain1 );
+ 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;
+ 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
+ idChainI = id;
+ 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");
+ 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; // 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) */
- atom1 = id; // choline or amine
- atom2 = id+1; // phosphate
- atom3 = id+2; // glycerol
- fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _phos_, atom1, atom2,
atom3 );
- break;
- case GLYCEROL:
- atom1 = id; // phosphate
- atom2 = id+1; // glycerol
- atom3 = id+2; // ester_a
- fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2,
atom3 ); // phos-gly-estA
- atom3 = id+8; // ester_b
- fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2,
atom3 ); // phos-gly-estB
- 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) */
- 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;
- 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+1;
- atom3 = id+2;
- fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, 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) */
- 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;
- default:
- printf("Error in function WriteLammpsData - default case reached.
Exiting...\n"); exit(0);
- }
+ for ( zRep = 0; zRep < zReplicas; zRep++ ) {
+ DO_SITE {
+ id = n + zRep * ( nAtoms / zReplicas ); // 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) */
+ atom1 = id; // choline or amine
+ atom2 = id+1; // phosphate
+ atom3 = id+2; // glycerol
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _phos_, atom1, atom2,
atom3 );
+ break;
+ case GLYCEROL:
+ atom1 = id; // phosphate
+ atom2 = id+1; // glycerol
+ atom3 = id+2; // ester_a
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2, atom3
); // phos-gly-estA
+ atom3 = id+8; // ester_b
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2, atom3
); // phos-gly-estB
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, DIPOLE, atom2, atom1,
atom1 ); // gly-phos angle_dipole
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, DIPOLE, id+2, id+3, id+3
); // estA-tailA1 angle_dipole
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, DIPOLE, id+8, id+9, id+9
); // estB-tailB1 angle_dipole
+ 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) */
+ 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;
+ 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+1;
+ atom3 = id+2;
+ fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, 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) */
+ 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;
+ 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);}
}
@@ -226,12 +236,19 @@
// fprintf(fp, "pair_style hybrid lj/cut 12 lj/cut/coul/cut 12
dipole/cut 9\n\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, "bond_style harmonic\n");
+ fprintf(fp, "angle_style hybrid cosine/squared dipole\n");
fprintf(fp, "\n");
for ( n = 0; n < nTypes-1; n++ ) { // write masses
- fprintf(fp, "mass%5d%9.3f\n", n+1, mass[n]);
+ if (n==1) {
+ if (nLipids==nDOPEs)
+ fprintf(fp, "mass%5d%9.3f\n", n+1, massAmine);
+ else if (nLipids==nDOPCsDSPCs)
+ fprintf(fp, "mass%5d%9.3f\n", n+1, mass[n]);
+ else {printf("Error: mass assignment for mixed lipid bilayers to be
coded.\n"); exit(0);}
+ } else
+ fprintf(fp, "mass%5d%9.3f\n", n+1, mass[n]);
}
fprintf(fp, "\n");
@@ -240,23 +257,44 @@
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);
- }
- }
- }
+ if ( i==1 || j==1 ){ // resolve AMINE case
+ if ( i==1 && j==1 ) {
+ if (nLipids==nDOPEs)
fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n",2,2,eps[AMINE_TYPE][AMINE_TYPE]*J_IN_cal,
sig[AMINE_TYPE][AMINE_TYPE]*nm_IN_Angstrom);
+ else if (nLipids==nDOPCsDSPCs)
fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n",2,2,eps[i][j]*J_IN_cal,
sig[i][j]*nm_IN_Angstrom);
+ else {printf("Error: pair coeff assignment for mixed lipid bilayers
to be coded.\n"); exit(0);}
+ } else {
+ if (j==1) {
+ if (nLipids==nDOPEs) fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n",
i+1, 2, eps[i][AMINE_TYPE]*J_IN_cal, sig[i][AMINE_TYPE]*nm_IN_Angstrom);
+ else if (nLipids==nDOPCsDSPCs)
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);
+ else {printf("Error: pair coeff assignment for mixed lipid bilayers
to be coded.\n"); exit(0);}
+ }
+ if (i==1) {
+ if (nLipids==nDOPEs) fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n",
i+1, j+1, eps[j][AMINE_TYPE]*J_IN_cal, sig[j][AMINE_TYPE]*nm_IN_Angstrom);
+ else if (nLipids==nDOPCsDSPCs)
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);
+ else {printf("Error: pair coeff assignment for mixed lipid bilayers
to be coded.\n"); exit(0);}
+ }
+ }
+ }
+ else 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");
if (nSites!=nWaters) { // to prevent segmentation fault when simulating
pure water systems
// 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 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
+ if (nLipids==nDOPCsDSPCs)
+ fprintf(fp, "bond_coeff%5d %7.3f%7.3f\n", X_PHOS,
springRigidity_LAMMPS_units/2,
0.9*sig[CHOLINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
+ else if (nLipids==nDOPEs)
+ fprintf(fp, "bond_coeff%5d %7.3f%7.3f\n", X_PHOS,
springRigidity_LAMMPS_units/2,
0.9*sig[AMINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
+ else {printf("Error: bond assignment for mixed lipid bilayers to be
coded.\n"); exit(0);}
+ 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", 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); */
@@ -268,10 +306,11 @@
// 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", _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 cosine/squared%8.3f%8.2f\n", _phos_,
angleRigidity_LAMMPS_units/2, refAngle_CholPhosGly*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d cosine/squared%8.3f%8.2f\n", _gly_,
angleRigidity_LAMMPS_units/2, refAngle_PhosGlyEst*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d cosine/squared%8.3f%8.2f\n", _piSat_,
angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d cosine/squared%8.3f%8.2f\n", _cisUnsat_,
angleRigidity_LAMMPS_units/2, refAngle_cisUnsat*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d dipole %8.3f%8.2f\n", DIPOLE,
orientRigidity_LAMMPS_units/2, 0.); // '0' is irrelevant at the moment, as
gamma_reference = 0
fprintf(fp, "\n");
}
}
=======================================
--- /trunk/src/main.c Wed Jul 20 07:59:24 2011
+++ /trunk/src/main.c Wed Sep 28 10:00:59 2011
@@ -185,8 +185,6 @@
real *rLipid; /* vector containing lipids Centre Of Geometry position */
extern VecI cells;
-extern int insertSolute; // if positive, activates solute insertion. The
actual value indicates how many solute molecules are inserted
-// in z-constraint context: logical input switch: if set to 1 a solute is
created and located at "zConstraintDist"
extern real rCutWatWat, rCutSolute, rCutSoluteElse, rCutLipLip;
extern const int applyBarostat, flexBox;
extern const int diffusion, stepDiffuse, latDiff, stepLatDiff,
writeLipLatMotion;
@@ -243,7 +241,6 @@
}
ConvertInputParametersToBrahmsUnits( fp );
- if ( zConstraint && insertSolute ) nSites++; // add one site - the
constrained solute
boxMatrix.u[0]=region.x; boxMatrix.u[4]=region.y;
boxMatrix.u[8]=region.z;
boxMatrix.u[1]=boxMatrix.u[2]=boxMatrix.u[3]=boxMatrix.u[5]=boxMatrix.u[6]=boxMatrix.u[7]=0.;
=======================================
--- /trunk/src/output.c Thu Aug 18 10:23:37 2011
+++ /trunk/src/output.c Wed Sep 28 10:00:59 2011
@@ -34,7 +34,7 @@
void PrintBrahmsVersion()
{
printf("\n**********************\n");
- printf( "* BRAHMS - 18Aug2011 *\n");
+ printf( "* BRAHMS - 27Sep2011 *\n");
printf( "**********************\n\n");
fflush(stdout);
}
=======================================
--- /trunk/src/startUp.c Tue Jun 28 02:10:55 2011
+++ /trunk/src/startUp.c Wed Sep 28 10:00:59 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.
@@ -32,7 +32,7 @@
extern real *rrDiffuseAv, *ddDiffuseAv, *latSqrDiffAv;
extern const int nValDiffuse, nBuffDiffuse, stepDiffuse, stepWatRdf,
limitWatRdf;
extern LatDiffBuff *latDiffBuff;
-extern const int nValLatDiff, nBuffLatDiff, stepLatDiff, limitLatDiffAv,
zConstraint, insertSolute;
+extern const int nValLatDiff, nBuffLatDiff, stepLatDiff, limitLatDiffAv,
zConstraint;
extern const int doCheckpoint, newRun, limitDiffuseAv, writeAreaVol,
flexBox;
extern real *histWatRdf, *totalHistEpp, *histLppNormalComponent,
*histLppLateralComponent;
extern real *histWpp, *soluteHistEpp, *watHistEpp, *headgroupHistEpp,
*glycerolHistEpp, *esterHistEpp;
@@ -759,15 +759,6 @@
if ( zConstraint || nSolutes ) FindSoluteNAtoms();
AllocArrays( fp );
- /* if ( newRun && loadStructure ) { */
- /* LoadStructureFile(); */
- /* if ( centerInputStruct ) CenterStructCoords(); */
- /* } */
- /* if ( zConstraint && nSolutes && !insertSolute ) { // if nSolutes ==
0, we are dealing with an SSD solute, so we don't have to
loadSoluteParameters */
- if ( zConstraint || nSolutes ) { // if nSolutes == 0, we are dealing
with an SSD solute, so we don't have to loadSoluteParameters
- LoadSoluteParameters( insertSolute, nAtoms );
- ConvertSoluteUnits();
- }
LoadTopology( fp );
InitializeMeasurements( fp );
if ( newRun ) {
@@ -808,8 +799,6 @@
xyArea.val = region.x * region.y;
ComputeTotalMass();
if ( resetTime ) { stepCount = 0; timeNow = 0.; }
- if ( zConstraint && insertSolute ) InsertSoluteForZConstraint();
- if ( nSolutes && insertSolute && !zConstraint )
InsertMultipleSolutes();
CheckStructureLoaded( nSites, nWaters );
printf( ">>>>>>>>>>>>>>>> CONTINUED RUN
<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" );
printf( "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n"
);
=======================================
--- /trunk/src/stepMain.c Tue Jan 18 06:19:42 2011
+++ /trunk/src/stepMain.c Wed Sep 28 10:00:59 2011
@@ -23,7 +23,7 @@
extern real timeNow;
extern Prop totEnergy, potEnergy, zConstraintForce;
extern const int applyThermostat, applyBarostat, flexBox,
removeSystemTranslation, removeMonolayersTranslation;
-extern const int nSites, nWaters, nLipids, nAtoms, nSolutes, DOFs,
zConstraint, insertSolute, stepNebr;
+extern const int nSites, nWaters, nLipids, nAtoms, nSolutes, DOFs,
zConstraint, stepNebr;
extern const int stepPdb, reCenterBilayer, stepLimit;
extern int stepCount;
extern RMat boxMatrix;
@@ -377,7 +377,6 @@
MZero( systemVirial.u ); /* initialise total system virial for a new
estimation at current step */
if ( nLipids ) StoreWaterCoorsPrevStep(); // for transbilayer water
permeability measurements
- if ( insertSolute && nSolutes ) GrowSolute(); /* increase solute's size
step by step */
LeapfrogStep( 1 );