Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Menu
Open sidebar
Development
openfoam
Commits
b33f76e5
Commit
b33f76e5
authored
Jun 23, 2008
by
Andrew Heather
Committed by
graham
Jun 23, 2008
Browse files
updating molecularDynamics functionality with Graham's latest changes
parent
12e0fc01
Changes
22
Hide whitespace changes
Inline
Side-by-side
applications/utilities/preProcessing/molConfig/Make/files
0 → 100755
View file @
b33f76e5
latticeStructures = latticeStructures
velocityDistributions = velocityDistributions
createMolecules.C
molConfig.C
genMolConfig.C
EXE = $(FOAM_APPBIN)/molConfig
applications/utilities/preProcessing/molConfig/correctVelocities.H
0 → 100644
View file @
b33f76e5
for
(
molN
=
totalMols
;
molN
<
totalMols
+
totalZoneMols
;
molN
++
)
{
// Remove bulk momentum introduced by random numbers and add
// desired bulk velocity
// For systems with molecules of significantly differing masses, this may
// need to be an iterative process or employ a better algorithm for
// removing an appropriate share of the excess momentum from each molecule.
initialVelocities
(
molN
)
+=
bulkVelocity
-
momentumSum
/
totalZoneMols
/
mass
;
}
// momentumSum = vector::zero;
//
// for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
// {
// momentumSum += mass*initialVelocities(molN);
// }
//
// Info << "Check momentum adjustment: " << momentumSum << endl;
applications/utilities/preProcessing/molConfig/createMolecules.C
0 → 100644
View file @
b33f76e5
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2005 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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 2 of the License, or (at your
option) any later version.
OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*----------------------------------------------------------------------------*/
#include
"molConfig.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void
Foam
::
molConfig
::
createMolecules
()
{
Info
<<
nl
<<
"Creating molecules from zone specifications
\n
"
<<
endl
;
DynamicList
<
vector
>
initialPositions
(
0
);
DynamicList
<
label
>
initialIds
(
0
);
DynamicList
<
scalar
>
initialMasses
(
0
);
DynamicList
<
label
>
initialCelli
(
0
);
DynamicList
<
vector
>
initialVelocities
(
0
);
DynamicList
<
vector
>
initialAccelerations
(
0
);
DynamicList
<
label
>
initialTethered
(
0
);
DynamicList
<
vector
>
initialTetherPositions
(
0
);
label
totalMols
=
0
;
label
idAssign
;
Random
rand
(
clock
::
getTime
());
// * * * * * * * * Building the IdList * * * * * * * * * //
//Notes: - each processor will have an identical idList_.
// - The order of id's inside the idList_ depends on the order
// of subDicts inside the molConigDict.
Info
<<
"Building the idList: "
;
forAll
(
molConfigDescription_
.
toc
(),
cZs
)
{
word
subDictName
(
molConfigDescription_
.
toc
()[
cZs
]);
word
iD
(
molConfigDescription_
.
subDict
(
subDictName
).
lookup
(
"id"
));
if
(
findIndex
(
idList_
,
iD
)
==
-
1
)
{
idList_
.
append
(
iD
);
}
}
forAll
(
idList_
,
i
)
{
Info
<<
" "
<<
idList_
[
i
];
}
Info
<<
nl
<<
endl
;
// * * * * * * * * Filling the Mesh * * * * * * * * * //
const
cellZoneMesh
&
cellZoneI
=
mesh_
.
cellZones
();
if
(
cellZoneI
.
size
())
{
Info
<<
"Filling the zones with molecules."
<<
nl
<<
endl
;
}
else
{
FatalErrorIn
(
"void createMolecules()
\n
"
)
<<
"No cellZones found in mesh description."
<<
abort
(
FatalError
);
}
forAll
(
cellZoneI
,
cZ
)
{
if
(
cellZoneI
[
cZ
].
size
())
{
if
(
!
molConfigDescription_
.
found
(
cellZoneI
[
cZ
].
name
()))
{
Info
<<
"Zone specification subDictionary: "
<<
cellZoneI
[
cZ
].
name
()
<<
" not found."
<<
endl
;
}
else
{
label
n
=
0
;
label
totalZoneMols
=
0
;
label
molsPlacedThisIteration
;
# include "readZoneSubDict.H"
idAssign
=
findIndex
(
idList_
,
id
);
# include "startingPoint.H"
// Continue trying to place molecules as long as at
// least one molecule is placed in each iteration.
// The "|| totalZoneMols == 0" condition means that the
// algorithm will continue if the origin is outside the
// zone - it will cause an infinite loop if no molecules
// are ever placed by the algorithm.
if
(
latticeStructure
!=
"empty"
)
{
while
(
molsPlacedThisIteration
!=
0
||
totalZoneMols
==
0
)
{
molsPlacedThisIteration
=
0
;
bool
partOfLayerInBounds
=
false
;
# include "createPositions.H"
if
(
totalZoneMols
==
0
&&
!
partOfLayerInBounds
)
{
WarningIn
(
"molConfig::createMolecules()"
)
<<
"A whole layer of unit cells was placed "
<<
"outside the bounds of the mesh, but no "
<<
"molecules have been placed in zone '"
<<
cellZoneI
[
cZ
].
name
()
<<
"'. This is likely to be because the zone "
<<
"has few cells ("
<<
cellZoneI
[
cZ
].
size
()
<<
" in this case) and no lattice position "
<<
"fell inside them. "
<<
"Aborting filling this zone."
<<
endl
;
break
;
}
totalZoneMols
+=
molsPlacedThisIteration
;
n
++
;
}
label
molN
;
for
(
molN
=
totalMols
;
molN
<
totalMols
+
totalZoneMols
;
molN
++
)
{
initialIds
.
append
(
idAssign
);
initialMasses
.
append
(
mass
);
initialAccelerations
.
append
(
vector
::
zero
);
if
(
tethered
)
{
initialTethered
.
append
(
1
);
initialTetherPositions
.
append
(
initialPositions
[
molN
]
);
}
else
{
initialTethered
.
append
(
0
);
initialTetherPositions
.
append
(
vector
::
zero
);
}
}
# include "createVelocities.H"
# include "correctVelocities.H"
}
totalMols
+=
totalZoneMols
;
}
}
}
idList_
.
shrink
();
positions_
=
initialPositions
;
positions_
.
setSize
(
initialPositions
.
size
());
id_
=
initialIds
;
id_
.
setSize
(
initialIds
.
size
());
mass_
=
initialMasses
;
mass_
.
setSize
(
initialMasses
.
size
());
cells_
=
initialCelli
;
cells_
.
setSize
(
initialCelli
.
size
());
U_
=
initialVelocities
;
U_
.
setSize
(
initialVelocities
.
size
());
A_
=
initialAccelerations
;
A_
.
setSize
(
initialAccelerations
.
size
());
tethered_
=
initialTethered
;
tethered_
.
setSize
(
initialTethered
.
size
());
tetherPositions_
=
initialTetherPositions
;
tetherPositions_
.
setSize
(
initialTetherPositions
.
size
());
nMol_
=
totalMols
;
}
// ************************************************************************* //
applications/utilities/preProcessing/molConfig/createPositions.H
0 → 100644
View file @
b33f76e5
vector
latticePosition
;
vector
globalPosition
;
if
(
latticeStructure
==
"SC"
)
{
# include "SC.H"
}
else
if
(
latticeStructure
==
"FCC"
)
{
# include "FCC.H"
}
else
if
(
latticeStructure
==
"BCC"
)
{
# include "BCC.H"
}
else
{
FatalErrorIn
(
"createPositions.H
\n
"
)
<<
"latticeStructure "
<<
latticeStructure
<<
" not supported."
<<
abort
(
FatalError
);
}
applications/utilities/preProcessing/molConfig/createVelocities.H
0 → 100644
View file @
b33f76e5
vector
velocity
;
vector
momentumSum
=
vector
::
zero
;
if
(
velocityDistribution
==
"uniform"
)
{
# include "uniform.H"
}
if
(
velocityDistribution
==
"maxwellian"
)
{
# include "maxwellian.H"
}
applications/utilities/preProcessing/molConfig/genMolConfig.C
0 → 100644
View file @
b33f76e5
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2005 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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 2 of the License, or (at your
option) any later version.
OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include
"molConfig.H"
#include
"fvCFD.H"
using
namespace
Foam
;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int
main
(
int
argc
,
char
*
argv
[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
Info
<<
nl
<<
"Reading molecular configuration description dictionary"
<<
endl
;
IOobject
molConfigDescriptionIOobject
(
"molConfigDict"
,
runTime
.
system
(),
runTime
,
IOobject
::
MUST_READ
,
IOobject
::
NO_WRITE
,
false
);
if
(
!
molConfigDescriptionIOobject
.
headerOk
())
{
FatalErrorIn
(
args
.
executable
())
<<
"Cannot find molConfig description file "
<<
nl
<<
args
.
caseName
()
/
runTime
.
system
()
/
"molConfig"
/
"molConfigDict"
<<
nl
<<
exit
(
FatalError
);
}
IOdictionary
molConfigDescription
(
molConfigDescriptionIOobject
);
// Create molCloud, registering object with mesh
Info
<<
nl
<<
"Creating molecular configuration"
<<
endl
;
molConfig
molecules
(
molConfigDescription
,
mesh
);
label
totalMolecules
=
molecules
.
nMol
();
if
(
Pstream
::
parRun
())
{
reduce
(
totalMolecules
,
sumOp
<
label
>
());
}
Info
<<
nl
<<
"Total number of molecules added: "
<<
totalMolecules
<<
nl
<<
endl
;
moleculeCloud
molCloud
(
mesh
,
molecules
.
nMol
(),
molecules
.
id
(),
molecules
.
mass
(),
molecules
.
positions
(),
molecules
.
cells
(),
molecules
.
U
(),
molecules
.
A
(),
molecules
.
tethered
(),
molecules
.
tetherPositions
()
);
IOdictionary
idListDict
(
IOobject
(
"idList"
,
mesh
.
time
().
constant
(),
mesh
,
IOobject
::
NO_READ
,
IOobject
::
AUTO_WRITE
)
);
idListDict
.
add
(
"idList"
,
molecules
.
molIdList
());
IOstream
::
defaultPrecision
(
12
);
Info
<<
nl
<<
"Writing molecular configuration"
<<
endl
;
if
(
!
mesh
.
write
())
{
FatalErrorIn
(
args
.
executable
())
<<
"Failed writing moleculeCloud."
<<
nl
<<
exit
(
FatalError
);
}
Info
<<
nl
<<
"ClockTime = "
<<
runTime
.
elapsedClockTime
()
<<
" s"
<<
nl
<<
endl
;
Info
<<
nl
<<
"End
\n
"
<<
endl
;
return
0
;
}
// ************************************************************************* //
applications/utilities/preProcessing/molConfig/latticeStructures/BCC.H
0 → 100644
View file @
b33f76e5
labelVector
iN
(
0
,
0
,
0
);
vector
gap
=
(
vector
::
one
)
*
pow
((
numberDensity
/
2
.
0
),
-
(
1
.
0
/
3
.
0
));
#include
"origin.H"
// Info<< "gap = " << gap << endl;
// Special treatment is required for the first position, i.e. iteration zero.
if
(
n
==
0
)
{
latticePosition
.
x
()
=
(
iN
.
x
()
*
gap
.
x
());
latticePosition
.
y
()
=
(
iN
.
y
()
*
gap
.
y
());
latticePosition
.
z
()
=
(
iN
.
z
()
*
gap
.
z
());
// Placing 2 molecules in each unit cell, using the algorithm from
// D. Rapaport, The Art of Molecular Dynamics Simulation, 2nd Ed, p68
for
(
label
iU
=
0
;
iU
<
2
;
iU
++
)
{
vector
unitCellLatticePosition
=
latticePosition
;
if
(
iU
==
1
)
{
unitCellLatticePosition
+=
0
.
5
*
gap
;
}
if
(
originSpecifies
==
"corner"
)
{
unitCellLatticePosition
-=
0
.
25
*
gap
;
}
// Info << nl << n << ", " << unitCellLatticePosition;
globalPosition
=
origin
+
transform
(
latticeToGlobal
,
unitCellLatticePosition
);
partOfLayerInBounds
=
mesh_
.
bounds
().
contains
(
globalPosition
);
if
(
findIndex
(
mesh_
.
cellZones
()[
cZ
],
mesh_
.
findCell
(
globalPosition
))
!=
-
1
)
{
molsPlacedThisIteration
++
;
initialPositions
.
append
(
globalPosition
);
initialCelli
.
append
(
mesh_
.
findCell
(
globalPosition
));
}
}
}
else
{
// Place top and bottom caps.
for
(
iN
.
z
()
=
-
n
;
iN
.
z
()
<=
n
;
iN
.
z
()
+=
2
*
n
)
{
for
(
iN
.
y
()
=
-
n
;
iN
.
y
()
<=
n
;
iN
.
y
()
++
)
{
for
(
iN
.
x
()
=
-
n
;
iN
.
x
()
<=
n
;
iN
.
x
()
++
)
{
latticePosition
.
x
()
=
(
iN
.
x
()
*
gap
.
x
());
latticePosition
.
y
()
=
(
iN
.
y
()
*
gap
.
y
());
latticePosition
.
z
()
=
(
iN
.
z
()
*
gap
.
z
());
for
(
label
iU
=
0
;
iU
<
2
;
iU
++
)
{
vector
unitCellLatticePosition
=
latticePosition
;
if
(
iU
==
1
)
{
unitCellLatticePosition
+=
0
.
5
*
gap
;
}
if
(
originSpecifies
==
"corner"
)
{
unitCellLatticePosition
-=
0
.
25
*
gap
;
}
// Info << nl << iN << ", " << unitCellLatticePosition;
globalPosition
=
origin
+
transform
(
latticeToGlobal
,
unitCellLatticePosition
);
partOfLayerInBounds
=
mesh_
.
bounds
().
contains
(
globalPosition
);
if
(
findIndex
(
mesh_
.
cellZones
()[
cZ
],
mesh_
.
findCell
(
globalPosition
)
)
!=
-
1
)
{
molsPlacedThisIteration
++
;
initialPositions
.
append
(
globalPosition
);
initialCelli
.
append
(
mesh_
.
findCell
(
globalPosition
));
}
}
}
}
}
// Placing sides
for
(
iN
.
z
()
=
-
(
n
-
1
);
iN
.
z
()
<=
(
n
-
1
);
iN
.
z
()
++
)
{
for
(
label
iR
=
0
;
iR
<=
2
*
n
-
1
;
iR
++
)
{
latticePosition
.
x
()
=
(
n
*
gap
.
x
());
latticePosition
.
y
()
=
((
-
n
+
(
iR
+
1
))
*
gap
.
y
());
latticePosition
.
z
()
=
(
iN
.
z
()
*
gap
.
z
());
for
(
label
iK
=
0
;
iK
<
4
;
iK
++
)
{
for
(
label
iU
=
0
;
iU
<
2
;
iU
++
)
{
vector
unitCellLatticePosition
=
latticePosition
;
if
(
iU
==
1
)
{
unitCellLatticePosition
+=
0
.
5
*
gap
;
}
if
(
originSpecifies
==
"corner"
)