[brahms-md] r48 committed - Added function 'WriteLammpsData' in 'output.c', and corresponding func...
brah...@googlecode.com
Author: orsimario
Date: Tue Jun 28 10:57:32 2011
Log: Added function 'WriteLammpsData' in 'output.c', and corresponding
function calls in 'main.c'. Changes to 'topology.c': added conversion of
timestep in reduced units for logging purposes, added auxiliary arrays to
store charges and dipoles, removed obsolete SSD-related stuff.
http://code.google.com/p/brahms-md/source/detail?r=48
Modified:
/trunk/src/dataStructs.h
/trunk/src/main.c
/trunk/src/output.c
/trunk/src/topology.c
=======================================
--- /trunk/src/dataStructs.h Thu Jan 6 06:46:39 2011
+++ /trunk/src/dataStructs.h Tue Jun 28 10:57:32 2011
@@ -102,7 +102,7 @@
typedef struct {
SiteType type; /* site type - see corresponding typedef */
MechType mechType; /* mechanical types - see corresponding typedef */
- int inLipid; // lipid identifier: -1 if site don't belong to any lipid
(eg water), 1 if it belongs to lipid-1, N if it belongs to lipid-N...
+ int inLipid; // lipid identifier: 0 if site don't belong to any lipid
(eg water), 1 if it belongs to lipid-1, N if it belongs to lipid-N...
LipidSpecies lipidSpecies; // lipid species (e.g., DOPC) see
corresponding typedef
LipidUnit lipidUnit; /* lipid unit (e.g., TAIL_A1) see corresponding
typedef */
LipidUnit skipNonbondedInteraction1; /* "first" unit excluded from
nonbonded interaction with current unit */
=======================================
--- /trunk/src/main.c Wed Feb 9 07:48:13 2011
+++ /trunk/src/main.c Tue Jun 28 10:57:32 2011
@@ -217,15 +217,15 @@
void PutCheckpoint();
void PutConfig();
void SetLipidTopologies();
-void SetSSDTopology();
void SetupJob( FILE* );
void SingleStep();
-void WriteAreaVolume();
void PrintNameLists( FILE* );
void PrintPdb( FILE* );
void GetNameLists( char** );
void SetupInterrupt();
void UnwrapLipLatCom();
+void WriteAreaVolume();
+void WriteLammpsData( FILE* );
/***********************************************************************
* SetParams -- sets various parameters; adapted from [Rapaport, p.31] *
@@ -235,7 +235,6 @@
const real WATER_NUMBER_DENSITY = 33.46; /* [molecules/nm^3] -
experimental number density of H20 */
SetLipidTopologies();
- SetSSDTopology();
if ( !nLipids && ( nWaters>255 ) ) { // build "pure" system of water
from scratch
VSCopy( region, 1. / pow( WATER_NUMBER_DENSITY / 4., 1./3. ),
initUcell );
@@ -282,8 +281,13 @@
SetupJob( fPtr );
fclose( fPtr );
+ fPtr = fopen( "data.initial", "w" );
+ WriteLammpsData( fPtr );
+ fclose( fPtr );
+
SetupInterrupt();
moreCycles = 1;
+
// PutConfig();
while ( moreCycles ) {
SingleStep();
@@ -331,6 +335,11 @@
if ( stepCount >= stepLimit ) moreCycles = 0;
} /* end while ( moreCycles ) */
if ( doCheckpoint ) PutCheckpoint(); /* put final checkpoint */
+
+ fPtr = fopen( "data.final", "w" );
+ WriteLammpsData( fPtr );
+ fclose( fPtr );
+
return 0; /* make the compiler happy */
}
=======================================
--- /trunk/src/output.c Tue Jun 28 02:10:55 2011
+++ /trunk/src/output.c Tue Jun 28 10:57:32 2011
@@ -19,11 +19,11 @@
extern const Atom *atom;
extern const LipidStruct dopcStruct;
extern const VecRTwo *xyLipidComTrue, xySoluteComTrue;
-extern const real timeNow, sigWat, sigChol, sigPhos, sigTail;
+extern const real timeNow, sigWat, watDipole, *mass, *charge, *dipole,
*sigma;
extern const Prop totEnergy, potEnergy, systemTemperature, lipTemperature,
watTemperature, solTemperature, pressure, dens;
extern const Prop boxVolume, xyArea, lipVolume, lipArea, surfaceTension,
headGroupDipole, headGroupTilt;
extern const Prop ordParTail, ordParTop, ordParMidUnsat, ordParMidSat,
ordParEnd, zConstraintForce;
-extern const int nSites, nWats, nLipids, nSolutes, stepCount, nAtoms,
zConstraint;
+extern const int nSites, nWats, nLipids, nSolutes, stepCount, nAtoms,
zConstraint, nTypes;
extern FILE *fPtrZCF, *fPtrAV;
/* external function declarations: */
@@ -45,11 +45,6 @@
printf( "% 8.3f % 8.3f % 8.3f\n", M.u[1], M.u[4], M.u[7] );
printf( "% 8.3f % 8.3f % 8.3f\n\n", M.u[2], M.u[5], M.u[8] );
}
-
-void WriteAreaVolume()
-{
- fprintf( fPtrAV, "%.3f %.2f\n", xyArea.val, boxVolume.val ); /* [nm^2],
[nm^3] */
-}
void PrintHeader( FILE *fp )
{
@@ -60,7 +55,7 @@
if ( nSolutes ) fprintf( fp, "%10s", "T-sol/C" );
fprintf( fp, "A_lip/(A^2) V_lip/(nm^3) ST/(dyn/cm) dipHG/D
tiltHG/deg" );
// fprintf(
fp, "%7s %7s %7s %8s %6s", "P2Top", "P2Mid", "P2Bot", "P2VBot", "Smol" );
-/* fprintf( fp, " %4s", "Smol" ); */
+ /* fprintf( fp, " %4s", "Smol" ); */
//fprintf( "%11s %10s ", "GlyTiltUp", "GlyTiltLo" );
fprintf( fp, "\n" );
} else { // if waters only
@@ -232,9 +227,9 @@
* causing a visualisation artefact. This bug does not affect the
simulation, but it would be good to fix it.
************************************************************************************************************/
static void UnwrapSitePositionWithRespectToGlycerol( VecR *rSiteWrapped,
/* out - wrapped position */
- int n, /* site number */
- int nSitesToGlycerol ) /* number of sites separating site "n" to
the
- glycerol site in the same (lipid) molecule */
+ int n, /* site number */
+ int nSitesToGlycerol ) /* number of sites separating site "n"
to the
+ glycerol site in the same (lipid) molecule */
{
VecR interSiteDistance;
VSub( interSiteDistance, siteCoords[ n ], siteCoords[ n-nSitesToGlycerol
] );
@@ -263,7 +258,8 @@
fprintf ( fPtr,"MODEL\n" );
// Adding unit cell information to display periodic images with VMD -
using PDB "CRYST1" record format (to be modified for non-tetragonal cells!):
- fprintf( fPtr, "CRYST1%9.3f%9.3f%9.3f 90.00 90.00 90.00 P 1 1 1
1\n", region.x*LENGTH_Angstrom, region.y*LENGTH_Angstrom,
region.z*LENGTH_Angstrom );
+ fprintf( fPtr, "CRYST1%9.3f%9.3f%9.3f 90.00 90.00 90.00 P 1 1 1
1\n",
+ region.x*LENGTH_Angstrom, region.y*LENGTH_Angstrom,
region.z*LENGTH_Angstrom );
DO_SITE {
rMat = siteRotMat[ n ];
VSCopy( rSite, LENGTH_Angstrom, siteCoords[n] );
@@ -288,7 +284,7 @@
fprintf( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "PO4", n+1, rSiteUnwrapped.x, rSiteUnwrapped.y,
rSiteUnwrapped.z );
break;
case GLYCEROL:
-/* fprintf( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "C", n+1, rSite.x, rSite.y, rSite.z ); */
+ /* fprintf( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "C", n+1, rSite.x, rSite.y, rSite.z ); */
VSAdd( tip, rSite, dipoleLength, siteOrientVec[ n ] );
fprintf ( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "C", n+1, tip.x, tip.y, tip.z );
VSAdd( tail, rSite, -dipoleLength, siteOrientVec[ n ] );
@@ -296,7 +292,7 @@
break;
case ESTER_A:
UnwrapSitePositionWithRespectToGlycerol( &rSiteUnwrapped, n, 1 );
- /* fprintf( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "O", n+1, rSiteUnwrapped.x, rSiteUnwrapped.y, rSiteUnwrapped.z
); */
+ /* fprintf( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "O", n+1, rSiteUnwrapped.x, rSiteUnwrapped.y, rSiteUnwrapped.z
); */
VSAdd( tip, rSiteUnwrapped, dipoleLength, siteOrientVec[ n ] );
fprintf ( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "O", n+1, tip.x, tip.y, tip.z );
VSAdd( tail,rSiteUnwrapped, -dipoleLength, siteOrientVec[ n ] );
@@ -304,7 +300,7 @@
break;
case ESTER_B: /* ester sites represented as 2-site molecules, to
highlight the dipole: "H" for the "+" end and "O" for the "-" */
UnwrapSitePositionWithRespectToGlycerol( &rSiteUnwrapped, n, 7 );
- /* fprintf( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "O", n+1, rSiteUnwrapped.x, rSiteUnwrapped.y, rSiteUnwrapped.z
); */
+ /* fprintf( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "O", n+1, rSiteUnwrapped.x, rSiteUnwrapped.y, rSiteUnwrapped.z
); */
VSAdd( tip, rSiteUnwrapped, dipoleLength, siteOrientVec[ n ] );
fprintf ( fPtr, "ATOM %6d %4s UNK %5d %11.3f%8.3f%8.3f\n",
++atomCount, "O", n+1, tip.x, tip.y, tip.z );
VSAdd( tail,rSiteUnwrapped, -dipoleLength, siteOrientVec[ n ] );
@@ -368,3 +364,49 @@
PrintPdb( fPtr );
fclose( fPtr );
}
+
+void WriteAreaVolume()
+{
+ fprintf( fPtrAV, "%.3f %.2f\n", xyArea.val, boxVolume.val ); /* [nm^2],
[nm^3] */
+}
+
+/****************************************************************************************************************
+ * WriteLammpsData -- writes "data." file that can be used by the LAMMPS
program
+
****************************************************************************************************************/
+void WriteLammpsData( FILE *fp )
+{
+ int n;
+
+ fprintf(fp, "Brahms-generated input file for LAMMPS\n");
+ fprintf(fp, "\n");
+
+ fprintf(fp, "%9d atoms\n", nSites);
+ fprintf(fp, "%9d atom types\n", nTypes);
+ 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");
+ fprintf(fp, "\n");
+ DO_SITE {
+ // write atom-ID and atom-type:
+ fprintf(fp, "%7d%5d", n+1, site[n].type+1 );
+ // write coordinates (x y z):
+ fprintf(fp, "%13.5f%13.5f%13.5f", siteCoords[n].x*LENGTH_Angstrom,
siteCoords[n].y*LENGTH_Angstrom, siteCoords[n].z*LENGTH_Angstrom);
+ // write molecule-ID (this is 0 for water sites):
+ fprintf(fp, "%6d", site[n].inLipid );
+ // write charge (this is 0.00 for water sites):
+ fprintf(fp, "%7.2f", charge[site[n].type]);
+ // write dipole component ('dipole' style):
+ fprintf(fp, "%13.5f%13.5f%13.5f",
+ siteOrientVec[ n ].x * dipole[site[n].type] * DIPOLE_D * D_IN_eA,
+ siteOrientVec[ n ].y * dipole[site[n].type] * DIPOLE_D * D_IN_eA,
+ siteOrientVec[ n ].z * dipole[site[n].type] * DIPOLE_D * D_IN_eA );
+ // write attributes for 'sphere' style - diameter & atom density:
+ fprintf(fp, "%7.1f%7.1f", sigma[site[n].type]*LENGTH_Angstrom,
mass[site[n].type] / (4*PI/3 *
Cube((sigma[site[n].type]/2)*LENGTH_Angstrom)));
+ fprintf(fp, "\n");
+ }
+}
=======================================
--- /trunk/src/topology.c Thu Jan 6 06:46:39 2011
+++ /trunk/src/topology.c Tue Jun 28 10:57:32 2011
@@ -14,7 +14,7 @@
extern const VecR region;
extern real **sig, **eps, **sigSquared, **productElectroMoments, *mass,
*inertia;
extern real sigWat, epsWat;
-extern const real extTemperature;
+extern const real extTemperature, deltaT;
extern const real sigChol, sigAmine, sigPhos, sigGly, sigEst, sigTail;
extern const real epsChol, epsAmine, epsPhos, epsEst, epsGly, epsTail;
extern const real massChol, massAmine, massPhos, massGly, soluteMass,
massEst, massTail, massWat;
@@ -28,7 +28,7 @@
const real hBondWatEst, hBondWatPhos, hBondWatGly, hBondWatWat;
const real hBondAmineWat, hBondAminePhos, hBondAmineEst, hBondAmineGly;
LipidStruct dopcStruct, dopeStruct;
-real *electroMoment, *typeSigma, *typeEpsilon;
+real *electroMoment, *sigma, *epsilon, *charge, *dipole;
real **productElectroMoments; /* product of electrostatic moments, i.e.
charge*charge charge*dipole or dipole*dipole */
VecR solInertia; // principal moments of inertia of the solute (rigid)
molecule
@@ -36,32 +36,11 @@
dopcStruct.nBeads = 15;
dopeStruct.nBeads = 15;
}
-
-void SetSSDTopology()
-{
- FILE *fp;
- // coordinates as in Bratko, Blum, Luzar, JCP, 1985
- /* real zOxygen = 0.00654; // [nm] hence already in brahms units */
- /* real zHydrogen = 0.053; // [nm] */
- /* real yHydrogen = 0.075; // [nm] */
- fp = fopen( "topology.log", "w" );
- /* fprintf( fp, "water Masses ( H ,O )= ( %3.1f, %3.1f )\n",
watHydrogenMass, watOxygenMass ); */
- fprintf( fp, "water mass = %2.0f\n", massWat );
- fprintf(fp,"water inertia = %5.4f u * nm^2\n", inertiaWat );
-
- /* watHydrogenMass = 15; // [amu] hence already in brahms units */
- /* watOxygenMass = 20; */
- /* ssdInertia.x = watOxygenMass * Sqr( zOxygen ) + 2. *
watHydrogenMass * ( Sqr( yHydrogen ) + Sqr( zHydrogen ) ); */
- /* ssdInertia.y = watOxygenMass * Sqr( zOxygen ) + 2. *
watHydrogenMass * Sqr( zHydrogen ); */
- /* ssdInertia.z = 2. * watHydrogenMass * Sqr( yHydrogen ); */
-
- fclose( fp );
-}
void LoadTopology( FILE *logFilePtr )
{
int type, i, j, counter = 1;
- real rTemperature, rDensity, rDipole;
+ real rTemperature, rDensity, rDipole, rTimeStep;
FILE *fp;
fprintf(logFilePtr, "LoadTopology... "); fflush(logFilePtr);
@@ -73,46 +52,62 @@
eps = CreateRealMatrix( nTypes, nTypes );
productElectroMoments = CreateRealMatrix( nTypes, nTypes );
- AllocMem( typeSigma, nTypes, real ); /* allocating arrays */
- AllocMem( typeEpsilon, nTypes, real );
+ AllocMem( sigma, nTypes, real ); /* allocating arrays */
+ AllocMem( epsilon, nTypes, real );
AllocMem( electroMoment, nTypes, real );
+ AllocMem( charge, nTypes, real );
+ AllocMem( dipole, nTypes, real );
for ( type = 0; type < nTypes; type++ ) { /* fill in arrays */
switch ( type ) {
case WATER_TYPE: /* SSD water */
- typeSigma[ type ] = sigWat;
- typeEpsilon[ type ] = epsWat;
+ sigma[ type ] = sigWat;
+ epsilon[ type ] = epsWat;
electroMoment[ type ] = watDipole;
+ charge[ type ] = 0;
+ dipole[ type ] = watDipole;
break;
case CHOLINE_TYPE: /* choline */
- typeSigma[ type ] = sigChol;
- typeEpsilon[ type ] = epsChol;
+ sigma[ type ] = sigChol;
+ epsilon[ type ] = epsChol;
electroMoment[ type ] = cholCharge;
+ charge[ type ] = cholCharge;
+ dipole[ type ] = 0;
break;
case AMINE_TYPE: /* amine */
- typeSigma[ type ] = sigAmine;
- typeEpsilon[ type ] = epsAmine;
+ sigma[ type ] = sigAmine;
+ epsilon[ type ] = epsAmine;
electroMoment[ type ] = cholCharge; /* amine (in DOPE) has the same
charge as choline (in DOPC) */
+ charge[ type ] = cholCharge;
+ dipole[ type ] = 0;
break;
case PHOSPHATE_TYPE: /* phosphate */
- typeSigma[ type ] = sigPhos;
- typeEpsilon[ type ] = epsPhos;
+ sigma[ type ] = sigPhos;
+ epsilon[ type ] = epsPhos;
electroMoment[ type ] = phosCharge;
+ charge[ type ] = phosCharge;
+ dipole[ type ] = 0;
break;
case GLYCEROL_TYPE:
- typeSigma[ type ] = sigGly;
- typeEpsilon[ type ] = epsGly;
+ sigma[ type ] = sigGly;
+ epsilon[ type ] = epsGly;
electroMoment[ type ] = glyDipole;
+ charge[ type ] = 0;
+ dipole[ type ] = glyDipole;
break;
case ESTER_TYPE:
- typeSigma[ type ] = sigEst;
- typeEpsilon[ type ] = epsEst;
+ sigma[ type ] = sigEst;
+ epsilon[ type ] = epsEst;
electroMoment[ type ] = estDipole;
+ charge[ type ] = 0;
+ dipole[ type ] = estDipole;
break;
case TAIL_TYPE:
- typeSigma[ type ] = sigTail;
- typeEpsilon[ type ] = epsTail;
+ sigma[ type ] = sigTail;
+ epsilon[ type ] = epsTail;
electroMoment[ type ] = 0;
+ charge[ type ] = 0;
+ dipole[ type ] = 0;
break;
case SOLUTE_TYPE: // temporary fix to avoid compilation error
break;
@@ -123,9 +118,9 @@
/* compile the interaction matrices according to Lorentz-Berthelot
mixing rules - ref: AllenTildesley, p.21 */
for ( i = 0; i < nTypes; i++ ) {
for ( j = 0; j < nTypes; j++ ) {
- sig[ i ][ j ] = 0.5 * ( typeSigma[ i ] + typeSigma[ j ] );
+ sig[ i ][ j ] = 0.5 * ( sigma[ i ] + sigma[ j ] );
sigSquared[ i ][ j ] = Sqr( sig[ i ][ j ] );
- eps[ i ][ j ] = sqrt( typeEpsilon[ i ] * typeEpsilon[ j ] );
+ eps[ i ][ j ] = sqrt( epsilon[ i ] * epsilon[ j ] );
/* electrostatics pair-magnitudes */
productElectroMoments[ i ][ j ] = electroMoment[ i ] *
electroMoment[ j ];
fprintf( fp, "item %6d: sig[ %d ][ %d ] = %.4f, sigSquared[ %d ][ %d
] = %.4f, eps[ %d ][ %d ] = %.4f \n",
@@ -138,6 +133,7 @@
if ( nWaters ) {
mass[ WATER_TYPE ] = massWat;
inertia[ WATER_TYPE ] = inertiaWat; // from input
+ fprintf( fp, "mass[ WATER_TYPE ] = %2.0f u\n", mass[ WATER_TYPE ] );
fprintf( fp, "inertia[ WATER_TYPE ] = %.4f u * nm^2\n", inertia[
WATER_TYPE ] );
/* WATER <--> WATER hydrogen-bonding */
eps[ WATER_TYPE ][ WATER_TYPE ] *= hBondWatWat;
@@ -153,8 +149,8 @@
mass[ ESTER_TYPE ] = massEst;
mass[ TAIL_TYPE ] = massTail;
/* computing principal moments of inertia */
-/* inertia[ GLYCEROL_TYPE ] = massGly * Sqr(sigGly) / 10; */
-/* inertia[ ESTER_TYPE ] = massEst * Sqr(sigEst) / 10; */
+ /* inertia[ GLYCEROL_TYPE ] = massGly * Sqr(sigGly) / 10; */
+ /* inertia[ ESTER_TYPE ] = massEst * Sqr(sigEst) / 10; */
inertia[ GLYCEROL_TYPE ] = inertiaGly;
inertia[ ESTER_TYPE ] = inertiaEst;
fprintf( fp, "inertia[ GLYCEROL_TYPE ] = %.4f u * nm^2\n", inertia[
GLYCEROL_TYPE ] );
@@ -207,14 +203,17 @@
}
fprintf( fp, "\nReduced values for water (Stockmayer fluid):\n" );
+ // From Appendix B of 'Computer simulation of liquids', Allen &
Tildesley, 1987
rTemperature = kB_IN_BRAHMS_UNITS * extTemperature /
(hBondWatWat*epsWat);
+ rDensity = (nWaters / (region.x*region.y*region.z)) * Cube( sigWat );
+ rDipole = watDipole / sqrt( 4 * PI * EPSILON0_IN_BRAHMS_UNITS *
(hBondWatWat*epsWat) * Cube( sigWat ) );
+ rTimeStep = sqrt( epsWat / ( mass[ WATER_TYPE ] * Sqr( sigWat ) ) ) *
deltaT;
fprintf( fp, "Temperature* = %.2f\n", rTemperature );
-
- rDensity = (nWaters / (region.x*region.y*region.z)) * Cube( sigWat );
fprintf( fp, "Density* = %.2f\n", rDensity );
-
- rDipole = watDipole / sqrt( 4 * PI * EPSILON0_IN_BRAHMS_UNITS *
(hBondWatWat*epsWat) * Cube( sigWat ) );
- fprintf( fp, "Dipole* = %.2f\n\n", rDipole );
+ fprintf( fp, "Dipole* = %.2f\n", rDipole );
+ fprintf( fp, "TimeStep* = %.4f\n\n", rTimeStep );
+
+
fclose( fp );
fprintf(logFilePtr, "done\n"); fflush(logFilePtr);