remove Fortran library based MESONT styles and the library itself
This commit is contained in:
Axel Kohlmeyer
@ -47,7 +47,6 @@ This is the list of packages that may require additional steps.
|
||||
* :ref:`LEPTON <lepton>`
|
||||
* :ref:`MACHDYN <machdyn>`
|
||||
* :ref:`MDI <mdi>`
|
||||
* :ref:`MESONT <mesont>`
|
||||
* :ref:`ML-HDNNP <ml-hdnnp>`
|
||||
* :ref:`ML-IAP <mliap>`
|
||||
* :ref:`ML-PACE <ml-pace>`
|
||||
@ -1852,48 +1851,6 @@ MDI package
|
||||
|
||||
----------
|
||||
|
||||
.. _mesont:
|
||||
|
||||
MESONT package
|
||||
-------------------------
|
||||
|
||||
This package includes a library written in Fortran 90 in the
|
||||
``lib/mesont`` folder, so a working Fortran 90 compiler is required to
|
||||
compile it. Also, the files with the force field data for running the
|
||||
bundled examples are not included in the source distribution. Instead
|
||||
they will be downloaded the first time this package is installed.
|
||||
|
||||
.. tabs::
|
||||
|
||||
.. tab:: CMake build
|
||||
|
||||
No additional settings are needed besides ``-D PKG_MESONT=yes``
|
||||
|
||||
.. tab:: Traditional make
|
||||
|
||||
Before building LAMMPS, you must build the *mesont* library in
|
||||
``lib/mesont``\ . You can also do it in one step from the
|
||||
``lammps/src`` dir, using a command like these, which simply
|
||||
invokes the ``lib/mesont/Install.py`` script with the specified
|
||||
args:
|
||||
|
||||
.. code-block:: bash
|
||||
|
||||
make lib-mesont # print help message
|
||||
make lib-mesont args="-m gfortran" # build with GNU g++ compiler (settings as with "make serial")
|
||||
make lib-mesont args="-m ifort" # build with Intel icc compiler
|
||||
|
||||
The build should produce two files: ``lib/mesont/libmesont.a`` and
|
||||
``lib/mesont/Makefile.lammps``\ . The latter is copied from an
|
||||
existing ``Makefile.lammps.\*`` and has settings needed to build
|
||||
LAMMPS with the *mesont* library (though typically the settings
|
||||
contain only the Fortran runtime library). If necessary, you can
|
||||
edit/create a new ``lib/mesont/Makefile.machine`` file for your
|
||||
system, which should define an ``EXTRAMAKE`` variable to specify a
|
||||
corresponding ``Makefile.lammps.machine`` file.
|
||||
|
||||
----------
|
||||
|
||||
.. _molfile:
|
||||
|
||||
MOLFILE package
|
||||
|
||||
@ -49,7 +49,6 @@ packages:
|
||||
* :ref:`LEPTON <lepton>`
|
||||
* :ref:`MACHDYN <machdyn>`
|
||||
* :ref:`MDI <mdi>`
|
||||
* :ref:`MESONT <mesont>`
|
||||
* :ref:`ML-HDNNP <ml-hdnnp>`
|
||||
* :ref:`ML-IAP <mliap>`
|
||||
* :ref:`ML-PACE <ml-pace>`
|
||||
|
||||
@ -88,7 +88,6 @@ KOKKOS, o = OPENMP, t = OPT.
|
||||
* :doc:`ke/atom/eff <compute_ke_atom_eff>`
|
||||
* :doc:`ke/eff <compute_ke_eff>`
|
||||
* :doc:`ke/rigid <compute_ke_rigid>`
|
||||
* :doc:`mesont <compute_mesont>`
|
||||
* :doc:`mliap <compute_mliap>`
|
||||
* :doc:`momentum <compute_momentum>`
|
||||
* :doc:`msd <compute_msd>`
|
||||
|
||||
@ -204,7 +204,6 @@ OPT.
|
||||
* :doc:`meam/sw/spline <pair_meam_sw_spline>`
|
||||
* :doc:`mesocnt <pair_mesocnt>`
|
||||
* :doc:`mesocnt/viscous <pair_mesocnt>`
|
||||
* :doc:`mesont/tpm <pair_mesont_tpm>`
|
||||
* :doc:`mgpt <pair_mgpt>`
|
||||
* :doc:`mie/cut (g) <pair_mie>`
|
||||
* :doc:`mliap (k) <pair_mliap>`
|
||||
|
||||
@ -59,6 +59,17 @@ Minimize style *fire/old* has been removed. Its functionality can be
|
||||
reproduced with *fire* with specific options. Please see the
|
||||
:doc:`min_modify command <min_modify>` documentation for details.
|
||||
|
||||
Pair style mesont/tpm, compute style mesont, atom style mesont
|
||||
--------------------------------------------------------------
|
||||
|
||||
.. deprecated:: TBD
|
||||
|
||||
Pair style *mesont/tpm*, compute style *mesont*, and atom style
|
||||
*mesont* have been removed from the :ref:`MESONT package <PKG-MESONT>`.
|
||||
The same functionality is available through
|
||||
:doc:`pair style mesocnt <pair_mesocnt>`,
|
||||
:doc:`bond style mesocnt <bond_mesocnt>` and
|
||||
:doc:`angle style mesocnt <angle_mesocnt>`.
|
||||
|
||||
REAX package
|
||||
------------
|
||||
|
||||
@ -1621,25 +1621,23 @@ MESONT package
|
||||
**Contents:**
|
||||
|
||||
MESONT is a LAMMPS package for simulation of nanomechanics of nanotubes
|
||||
(NTs). The model is based on a coarse-grained representation of NTs as
|
||||
(NTs). The model is based on a coarse-grained representation of NTs as
|
||||
"flexible cylinders" consisting of a variable number of
|
||||
segments. Internal interactions within a NT and the van der Waals
|
||||
interaction between the tubes are described by a mesoscopic force field
|
||||
designed and parameterized based on the results of atomic-level
|
||||
molecular dynamics simulations. The description of the force field is
|
||||
provided in the papers listed below.
|
||||
provided in the papers listed in ``src/MESONT/README``.
|
||||
|
||||
This package contains two independent implementations of this model:
|
||||
:doc:`pair_style mesont/tpm <pair_mesont_tpm>` is the original
|
||||
implementation of the model based on a Fortran library in the
|
||||
``lib/mesont`` folder. The second implementation is provided by the
|
||||
mesocnt styles (:doc:`bond_style mesocnt <bond_mesocnt>`,
|
||||
:doc:`angle_style mesocnt <angle_mesocnt>` and :doc:`pair_style mesocnt
|
||||
<pair_mesocnt>`). The mesocnt implementation has the same features as
|
||||
the original implementation with the addition of friction, but is
|
||||
directly implemented in C++, interfaces more cleanly with general LAMMPS
|
||||
functionality, and is typically faster. It also does not require its own
|
||||
atom style and can be installed without any external libraries.
|
||||
This package used to have two independent implementations of this model:
|
||||
the original implementation using a Fortran library written by the
|
||||
developers of the model and a second implementation written in C++ by
|
||||
Philipp Kloza (U Cambridge). Since the C++ implementation offers the
|
||||
same features as the original implementation with the addition of
|
||||
friction, is typically faster, and easier to compile/install, the
|
||||
Fortran library based implementation has since been obsoleted and
|
||||
removed from the distribution. You have to download and compile
|
||||
an older version of LAMMPS if you want to use those.
|
||||
|
||||
**Download of potential files:**
|
||||
|
||||
@ -1648,12 +1646,14 @@ not included in the regular downloaded packages of LAMMPS or the git
|
||||
repositories. Instead, they will be automatically downloaded from a web
|
||||
server when the package is installed for the first time.
|
||||
|
||||
**Authors of the *mesont* styles:**
|
||||
**Authors of the obsoleted *mesont* styles:**
|
||||
|
||||
Maxim V. Shugaev (University of Virginia), Alexey N. Volkov (University
|
||||
of Alabama), Leonid V. Zhigilei (University of Virginia)
|
||||
|
||||
**Author of the *mesocnt* styles:**
|
||||
.. deprecated:: TBD
|
||||
|
||||
**Author of the C++ styles:**
|
||||
Philipp Kloza (U Cambridge)
|
||||
|
||||
.. versionadded:: 15Jun2020
|
||||
@ -1662,14 +1662,10 @@ Philipp Kloza (U Cambridge)
|
||||
|
||||
* src/MESONT: filenames -> commands
|
||||
* src/MESONT/README
|
||||
* :doc:`atom_style mesont <atom_style>`
|
||||
* :doc:`pair_style mesont/tpm <pair_mesont_tpm>`
|
||||
* :doc:`compute mesont <compute_mesont>`
|
||||
* :doc:`bond_style mesocnt <bond_mesocnt>`
|
||||
* :doc:`angle_style mesocnt <angle_mesocnt>`
|
||||
* :doc:`pair_style mesocnt <pair_mesocnt>`
|
||||
* examples/PACKAGES/mesont
|
||||
* tools/mesont
|
||||
|
||||
----------
|
||||
|
||||
|
||||
@ -275,9 +275,9 @@ whether an extra library is needed to build and use the package:
|
||||
- no
|
||||
* - :ref:`MESONT <PKG-MESONT>`
|
||||
- mesoscopic tubular potential model
|
||||
- pair styles :doc:`mesont/tpm <pair_mesont_tpm>`, :doc:`mesocnt <pair_mesocnt>`
|
||||
- pair styles :doc:`mesocnt <pair_mesocnt>`
|
||||
- PACKAGES/mesont
|
||||
- int
|
||||
- no
|
||||
* - :ref:`MGPT <PKG-MGPT>`
|
||||
- fast MGPT multi-ion potentials
|
||||
- :doc:`pair_style mgpt <pair_mgpt>`
|
||||
|
||||
@ -10,7 +10,7 @@ Syntax
|
||||
|
||||
atom_style style args
|
||||
|
||||
* style = *amoeba* or *angle* or *atomic* or *body* or *bond* or *charge* or *dielectric* or *dipole* or *dpd* or *edpd* or *electron* or *ellipsoid* or *full* or *line* or *mdpd* or *mesont* or *molecular* or *oxdna* or *peri* or *smd* or *sph* or *sphere* or *bpm/sphere* or *spin* or *tdpd* or *tri* or *template* or *wavepacket* or *hybrid*
|
||||
* style = *amoeba* or *angle* or *atomic* or *body* or *bond* or *charge* or *dielectric* or *dipole* or *dpd* or *edpd* or *electron* or *ellipsoid* or *full* or *line* or *mdpd* or *molecular* or *oxdna* or *peri* or *smd* or *sph* or *sphere* or *bpm/sphere* or *spin* or *tdpd* or *tri* or *template* or *wavepacket* or *hybrid*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
@ -109,8 +109,6 @@ quantities.
|
||||
+--------------+-----------------------------------------------------+--------------------------------------+
|
||||
| *mdpd* | density | mDPD particles |
|
||||
+--------------+-----------------------------------------------------+--------------------------------------+
|
||||
| *mesont* | mass, radius, length, buckling, connections, tube id| mesoscopic nanotubes |
|
||||
+--------------+-----------------------------------------------------+--------------------------------------+
|
||||
| *molecular* | bonds, angles, dihedrals, impropers | uncharged molecules |
|
||||
+--------------+-----------------------------------------------------+--------------------------------------+
|
||||
| *oxdna* | nucleotide polarity | coarse-grained DNA and RNA models |
|
||||
@ -382,8 +380,6 @@ The *sph* style is part of the SPH package for smoothed particle
|
||||
hydrodynamics (SPH). See `this PDF guide
|
||||
<PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in LAMMPS.
|
||||
|
||||
The *mesont* style is part of the MESONT package.
|
||||
|
||||
The *spin* style is part of the SPIN package.
|
||||
|
||||
The *wavepacket* style is part of the AWPMD package for the
|
||||
|
||||
@ -254,7 +254,6 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
|
||||
* :doc:`pair/local <compute_pair_local>` - distance/energy/force of each pairwise interaction
|
||||
* :doc:`pe <compute_pe>` - potential energy
|
||||
* :doc:`pe/atom <compute_pe_atom>` - potential energy for each atom
|
||||
* :doc:`mesont <compute_mesont>` - Nanotube bending,stretching, and intertube energies
|
||||
* :doc:`pe/mol/tally <compute_tally>` - potential energy between two groups of atoms separated into intermolecular and intramolecular components via the tally callback mechanism
|
||||
* :doc:`pe/tally <compute_tally>` - potential energy between two groups of atoms via the tally callback mechanism
|
||||
* :doc:`plasticity/atom <compute_plasticity_atom>` - Peridynamic plasticity for each atom
|
||||
|
||||
@ -1,59 +0,0 @@
|
||||
.. index:: compute mesont
|
||||
|
||||
compute mesont command
|
||||
======================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute ID group-ID mesont style
|
||||
|
||||
* ID, group-ID are documented in :doc:`compute <compute>` command
|
||||
* mesont = style name of the compute command
|
||||
* style = *estretch* or *ebend* or *etube*
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute 1 all mesont estretch
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
These computes define computations for the stretching (*estretch*), bending
|
||||
(*ebend*), and intertube (*etube*) per-node (atom) and total energies. The
|
||||
evaluated value is selected by the style parameter passed to the compute
|
||||
(*estretch*, *ebend*,or *etube*).
|
||||
|
||||
Output info
|
||||
"""""""""""
|
||||
|
||||
These computes calculate per-node (per-atom) vectors, which can be accessed by
|
||||
any command that uses per-atom values from a compute as input, and global
|
||||
scalars. See the :doc:`Howto output <Howto_output>` page for an overview of
|
||||
LAMMPS output options.
|
||||
|
||||
The computed values are provided in energy :doc:`units <units>`.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
These computes are part of the MESONT package. They are only enabled if
|
||||
LAMMPS is built with that package. See the :doc:`Build package <Build_package>`
|
||||
doc page for more info. In addition, :doc:`mesont pair_style <pair_style>`
|
||||
must be used.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`dump custom <dump>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
none
|
||||
|
||||
@ -90,11 +90,6 @@ Syntax
|
||||
de = change in thermal energy
|
||||
cv = heat capacity
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
MESONT package per-atom properties:
|
||||
buckling = buckling flag used in mesoscopic simulation of nanotubes
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
|
||||
@ -216,9 +216,8 @@ LAMMPS was built with that package. See the :doc:`Build package
|
||||
These pair styles require the :doc:`newton <newton>` setting to be
|
||||
"on" for pair interactions.
|
||||
|
||||
These pair styles require all 3 :doc:`special_bonds lj
|
||||
<special_bonds>` settings to be non-zero for proper neighbor list
|
||||
construction.
|
||||
These pair styles require all 3 :doc:`special_bonds lj <special_bonds>`
|
||||
settings to be non-zero for proper neighbor list construction.
|
||||
|
||||
Pair style *mesocnt/viscous* requires you to use the :doc:`comm_modify
|
||||
vel yes <comm_modify>` command so that velocities are stored by ghost
|
||||
@ -227,7 +226,9 @@ atoms.
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`pair_coeff <pair_coeff>`
|
||||
:doc:`pair_coeff <pair_coeff>`,
|
||||
:doc:`bond_style mesocnt <bond_mesocnt>`,
|
||||
:doc:`angle_style mesocnt <angle_mesocnt>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
@ -1,219 +0,0 @@
|
||||
.. index:: pair_style mesont/tpm
|
||||
|
||||
pair_style mesont/tpm command
|
||||
=============================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
pair_style mesont/tpm cut table_path BendingMode TPMType
|
||||
|
||||
* cut = the cutoff distance
|
||||
* table_path = the path to the potential table
|
||||
* BendingMode = the parameter defining the type of the bending potential for nanotubes: 0 - harmonic bending :ref:`(Srivastava) <Srivastava>`, 1 - anharmonic potential of bending and bending-buckling :ref:`(Zhigilei1) <Zhigilei1>`
|
||||
* TPMType = the parameter determining the type of the inter-tube interaction term: 0 - segment-segment approach, 1 - segment-chain approach :ref:`(Zhigilei2 <Zhigilei2>`, :ref:`Zhigilei3) <Zhigilei3>`
|
||||
|
||||
The segment-segment approach is approximately 5 times slower than segment-chain approximation.
|
||||
The parameter BendingMode also affects the calculation of the inter-tube interaction term when TPMType = 1. In this case, when BendingMode = 1, each continuous chain of segments is additionally replaced by a number of sub-chains divided by bending buckling kinks.
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
pair_style mesont/tpm 30.0 MESONT-TABTP_10_10.xrs 0 0
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
The tubular potential model (TPM) force field is designed for mesoscopic
|
||||
simulations of interacting flexible nanotubes. The force field is based on the
|
||||
mesoscopic computational model suggested in Ref. :ref:`(Srivastava) <Srivastava>`.
|
||||
In this model, each nanotube is represented by a chain of mesoscopic elements
|
||||
in the form of stretchable cylindrical segments, where each segment consists
|
||||
of multiple atoms. Each nanotube is divided into segments by a sequence of
|
||||
nodes placed on the nanotube centerline. This sequence of nodes determines the
|
||||
spatial position of the cylindrical segments and defines the configuration of
|
||||
the entire tube.
|
||||
|
||||
The potential force field that controls the dynamic behavior of a system of
|
||||
interacting nanotubes is given by the following equation defining the potential
|
||||
energy of the system:
|
||||
|
||||
.. math::
|
||||
|
||||
U = U_{str} + U_{bnd} + U_{vdW}
|
||||
|
||||
where :math:`U_{str}` is the harmonic potential describing the stretching of nanotube
|
||||
:ref:`(Srivastava) <Srivastava>`, :math:`U_{bnd}` is the potential for nanotube bending
|
||||
:ref:`(Srivastava) <Srivastava>` and bending-buckling :ref:`(Zhigilei1) <Zhigilei1>`, and
|
||||
:math:`U_{vdW}` is the potential describing van-der Waals interaction between nanotubes
|
||||
:ref:`(Zhigilei2 <Zhigilei2>`, :ref:`Zhigilei3) <Zhigilei3>`. The stretching energy, :math:`U_{str}` ,
|
||||
is given by the sum of stretching energies of individual nanotube segments.
|
||||
The bending energy, :math:`U_{bnd}` , is given by the sum of bending energies in all
|
||||
internal nanotube nodes. The tube-tube interaction energy, :math:`U_{vdW}` , is calculated
|
||||
based on the tubular potential method suggested in Ref. :ref:`(Zhigilei2) <Zhigilei2>`.
|
||||
The tubular potential method is briefly described below.
|
||||
|
||||
The interaction between two straight nanotubes of arbitrary length and
|
||||
orientation is described by the approximate tubular potential developed in
|
||||
:ref:`(Zhigilei3) <Zhigilei3>`. This potential approximates the results of direct
|
||||
integration of carbon-carbon interatomic potential over the surfaces of the
|
||||
interacting nanotubes, with the force sources homogeneously distributed over
|
||||
the nanotube surfaces. The input data for calculation of tubular potentials
|
||||
are partially tabulated. For single-walled CNTs of arbitrary chirality, the
|
||||
tabulated potential data can be generated in the form of ASCII files
|
||||
TPMSSTP.xrs and TPMA.xrs by the stand-alone code TMDPotGen included in the
|
||||
tool directory of LAMMPS release. The potential provided with LAMMPS release,
|
||||
MESONT-TABTP_10_10.xrs, is tabulated for (10,10) nanotubes.
|
||||
|
||||
Calculations of the interaction between curved or bent nanotubes are performed
|
||||
on either segment-segment or segment-chain basis. In the first case, activated
|
||||
when parameter TPMType is equal to 0, the tubular potential is calculated for
|
||||
each pair of interacting mesoscopic segments. In this case, however, small
|
||||
potential barriers for inter-tube sliding are introduced. While relatively
|
||||
small, these barriers are still larger than the ones that originate from the
|
||||
atomic-scale corrugation in atomistic modeling of inter-tube interaction. The
|
||||
latter are too weak to prevent room-temperature rearrangements of defect-free
|
||||
CNT, while the artificial mesoscopic barriers due to the segment-segment
|
||||
interaction can impede sliding of nanotubes with respect to each other and
|
||||
affect the kinetics of structural rearrangements in a system of nanotubes at
|
||||
moderate mesoscopic temperatures. In the second case, activated when parameter
|
||||
TPMType is equal to 1, the inter-tube interaction term is calculated based on
|
||||
the segment-chain approach. In this case, for each NT segment, the list of its
|
||||
neighboring segments is divided into short continuous chains of segments
|
||||
belonging to individual nanotubes. For each pair of a segment and a chain, the
|
||||
curved chain is approximated by a straight equivalent nanotube based on the
|
||||
weighted approach suggested in Ref. :ref:`(Zhigilei2) <Zhigilei2>`. Finally, the
|
||||
interaction between the segment and straight equivalent chain is calculated
|
||||
based on the tubular potential. In this case, and in the absence of bending
|
||||
buckling (i.e., when parameter BendingMode is equal to 0), the tubular
|
||||
potential method ensures the absence of corrugation of the effective inter-tube
|
||||
interaction potential for curved nanotubes and eliminates any barriers for the
|
||||
inter-tube sliding. As a result, the tubular potential method can describe the
|
||||
spontaneous self-assembly of nanotubes into continuous networks of bundles
|
||||
:ref:`(Zhigilei1 <Zhigilei1>`, :ref:`Zhigilei3) <Zhigilei3>`.
|
||||
|
||||
|
||||
----------
|
||||
|
||||
|
||||
The TMD force field has been used for generation of nanotube films, fibers,
|
||||
and vertically aligned forests of nanotubes. Mesoscopic dynamic simulations
|
||||
were used to prepare realistic structures of continuous networks of nanotube
|
||||
bundles and to study their structural and mechanical properties
|
||||
:ref:`(Zhigilei1 <Zhigilei1>`, :ref:`Zhigilei3 <Zhigilei3>`, :ref:`Zhigilei4 <Zhigilei4>`,
|
||||
:ref:`Zhigilei5 <Zhigilei5>`, :ref:`Zhigilei6) <Zhigilei6>`. With
|
||||
additional models for heat transfer, this force filed was also used to
|
||||
study the thermal transport properties of carbon nanotube films
|
||||
:ref:`(Zhigilei7 <Zhigilei7>`, :ref:`Zhigilei8 <Zhigilei8>`, :ref:`Zhigilei8) <Zhigilei8>`.
|
||||
The methods for modeling of
|
||||
the mechanical energy dissipation into heat (energy exchange between the
|
||||
dynamic degrees of freedom of the mesoscopic model and the energy of atomic
|
||||
vibrations that are not explicitly represented in the model)
|
||||
:ref:`(Zhigilei10) <Zhigilei10>` and mesoscopic description of covalent cross-links
|
||||
between nanotubes :ref:`(Banna) <Banna>` have also been developed but are not
|
||||
included in this first release of the LAMMPS implementation of the force field.
|
||||
Further details can be found in references provided below.
|
||||
|
||||
The MESONT package also provides TMDGen code designed to generate initial samples
|
||||
composed of straight and dispersed nanotubes of given chirality and length at a
|
||||
given material density, which is available in tools directory. In the generated
|
||||
samples, nanotubes are distributed with random positions and orientations. Both
|
||||
periodic and free boundary conditions are available along each axis of the
|
||||
system of coordinates. All parameters in the sample files generated with TMDGen
|
||||
are given in metal :doc:`units <units>`.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
|
||||
This pair style is a part of the MSEONT package, and it is only enabled if
|
||||
LAMMPS is built with that package. See the :doc:`Build package <Build_package>`
|
||||
doc page for more information.
|
||||
|
||||
This pair potential requires use of :doc:`mesont atomic style <atom_style>`.
|
||||
|
||||
This pair potential requires the :doc:`newton <newton>` setting to be "on" for
|
||||
pair interactions.
|
||||
|
||||
The cutoff distance should be set to be at least :math:`max\left[2L,\sqrt{L^2/2+(2R+T_{cut})^2}\right]` ,
|
||||
where L is the maximum segment length, R is the maximum tube radius, and
|
||||
:math:`T_{cut}` = 10.2 A is the maximum distance between the surfaces of interacting
|
||||
segments. Because of the use of extended chain concept at CNT ends, the recommended
|
||||
cutoff is 3L.
|
||||
|
||||
.. note::
|
||||
|
||||
Because of their size, *mesont* style potential files
|
||||
are not bundled with LAMMPS. When compiling LAMMPS from
|
||||
source code, the file ``TABTP_10_10.mesont`` should be downloaded
|
||||
transparently from `https://download.lammps.org/potentials/TABTP_10_10.mesont <https://download.lammps.org/potentials/TABTP_10_10.mesont>`_
|
||||
|
||||
The ``TABTP_10_10.mesont`` potential file is parameterized for metal :doc:`units <units>`.
|
||||
You can use the carbon nanotube mesoscopic force field with any LAMMPS units,
|
||||
but you would need to create your own potential files with coefficients listed in
|
||||
appropriate units, if your simulation does not use "metal" units.
|
||||
|
||||
The chirality parameters set during system generation must match the values
|
||||
specified during generation of the potential tables.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`pair_coeff <pair_coeff>`
|
||||
|
||||
----------
|
||||
|
||||
.. _Srivastava:
|
||||
|
||||
**(Srivastava)** Zhigilei, Wei, Srivastava, Phys. Rev. B 71, 165417 (2005).
|
||||
|
||||
.. _Zhigilei1:
|
||||
|
||||
**(Zhigilei1)** Volkov and Zhigilei, ACS Nano 4, 6187 (2010).
|
||||
|
||||
.. _Zhigilei2:
|
||||
|
||||
**(Zhigilei2)** Volkov, Simov, Zhigilei, ASME paper IMECE2008, 68021 (2008).
|
||||
|
||||
.. _Zhigilei3:
|
||||
|
||||
**(Zhigilei3)** Volkov, Zhigilei, J. Phys. Chem. C 114, 5513 (2010).
|
||||
|
||||
.. _Zhigilei4:
|
||||
|
||||
**(Zhigilei4)** Wittmaack, Banna, Volkov, Zhigilei, Carbon 130, 69 (2018).
|
||||
|
||||
.. _Zhigilei5:
|
||||
|
||||
**(Zhigilei5)** Wittmaack, Volkov, Zhigilei, Compos. Sci. Technol. 166, 66 (2018).
|
||||
|
||||
.. _Zhigilei6:
|
||||
|
||||
**(Zhigilei6)** Wittmaack, Volkov, Zhigilei, Carbon 143, 587 (2019).
|
||||
|
||||
.. _Zhigilei7:
|
||||
|
||||
**(Zhigilei7)** Volkov, Zhigilei, Phys. Rev. Lett. 104, 215902 (2010).
|
||||
|
||||
.. _Zhigilei8:
|
||||
|
||||
**(Zhigilei8)** Volkov, Shiga, Nicholson, Shiomi, Zhigilei, J. Appl. Phys. 111, 053501 (2012).
|
||||
|
||||
.. _Zhigilei9:
|
||||
|
||||
**(Zhigilei9)** Volkov, Zhigilei, Appl. Phys. Lett. 101, 043113 (2012).
|
||||
|
||||
.. _Zhigilei10:
|
||||
|
||||
**(Zhigilei10)** Jacobs, Nicholson, Zemer, Volkov, Zhigilei, Phys. Rev. B 86, 165414 (2012).
|
||||
|
||||
.. _Banna:
|
||||
|
||||
**(Banna)** Volkov, Banna, Comp. Mater. Sci. 176, 109410 (2020).
|
||||
|
||||
@ -283,7 +283,6 @@ accelerated styles exist.
|
||||
* :doc:`mesocnt <pair_mesocnt>` - mesoscopic vdW potential for (carbon) nanotubes
|
||||
* :doc:`mesocnt/viscous <pair_mesocnt>` - mesoscopic vdW potential for (carbon) nanotubes with friction
|
||||
* :doc:`mgpt <pair_mgpt>` - simplified model generalized pseudopotential theory (MGPT) potential
|
||||
* :doc:`mesont/tpm <pair_mesont_tpm>` - nanotubes mesoscopic force field
|
||||
* :doc:`mie/cut <pair_mie>` - Mie potential
|
||||
* :doc:`mm3/switch3/coulgauss/long <pair_lj_switch3_coulgauss_long>` - smoothed MM3 vdW potential with Gaussian electrostatics
|
||||
* :doc:`momb <pair_momb>` - Many-Body Metal-Organic (MOMB) force field
|
||||
|
||||
@ -710,8 +710,6 @@ of analysis.
|
||||
- atom-ID molecule-ID atom-type lineflag density x y z
|
||||
* - mdpd
|
||||
- atom-ID atom-type rho x y z
|
||||
* - mesont
|
||||
- atom-ID molecule-ID atom-type bond_nt mass mradius mlength buckling x y z
|
||||
* - molecular
|
||||
- atom-ID molecule-ID atom-type x y z
|
||||
* - peri
|
||||
@ -742,8 +740,6 @@ The per-atom values have these meanings and units, listed alphabetically:
|
||||
* atom-ID = integer ID of atom
|
||||
* atom-type = type of atom (1-Ntype, or type label)
|
||||
* bodyflag = 1 for body particles, 0 for point particles
|
||||
* bond_nt = bond NT factor for MESONT particles (?? units)
|
||||
* buckling = buckling factor for MESONT particles (?? units)
|
||||
* ccN = chemical concentration for tDPD particles for each species (mole/volume units)
|
||||
* cradius = contact radius for SMD particles (distance units)
|
||||
* cs_re,cs_im = real/imaginary parts of wave packet coefficients
|
||||
@ -760,9 +756,7 @@ The per-atom values have these meanings and units, listed alphabetically:
|
||||
* kradius = kernel radius for SMD particles (distance units)
|
||||
* lineflag = 1 for line segment particles, 0 for point or spherical particles
|
||||
* mass = mass of particle (mass units)
|
||||
* mlength = ?? length for MESONT particles (distance units)
|
||||
* molecule-ID = integer ID of molecule the atom belongs to
|
||||
* mradius = ?? radius for MESONT particles (distance units)
|
||||
* mux,muy,muz = components of dipole moment of atom (dipole units)
|
||||
* q = charge on atom (charge units)
|
||||
* rho = density (need units) for SPH particles
|
||||
|
||||
@ -4,16 +4,7 @@
|
||||
The files in this folder provide examples of using the CNT
|
||||
mesoscopic force field (MESONT).
|
||||
|
||||
Contributing author: Maxim Shugaev (UVA), mvs9t@virginia.edu
|
||||
|
||||
"bundle" is an example with a single bundle composed of 7 nanotubes
|
||||
using the mesont/tpm pair style
|
||||
|
||||
"film" is an example with a film composed of 396 200-nm-long
|
||||
nanotubes (79596 nodes) using the mesont/tpm pair style
|
||||
|
||||
|
||||
Contributing author: Philipp Kloza (U Cambridge), pak37@cam.ac.uk
|
||||
|
||||
"cnt" is an example showing CNT aerogel formation
|
||||
"film_mesocnt" is an example showing CNT aerogel formation
|
||||
using the mesocnt pair style
|
||||
|
||||
@ -1 +0,0 @@
|
||||
../../../potentials/TABTP_10_10.mesont
|
||||
@ -1,93 +0,0 @@
|
||||
|
||||
|
||||
77 atoms
|
||||
|
||||
1 atom types
|
||||
|
||||
-143.89 143.89 xlo xhi
|
||||
-143.89 143.89 ylo yhi
|
||||
0 220 zlo zhi
|
||||
|
||||
Masses
|
||||
|
||||
1 1.0
|
||||
|
||||
Atoms
|
||||
|
||||
1 1 1 11 2 5860.43 6.785 20 0 0 0 0 0 0 0
|
||||
2 1 1 1 3 5860.43 6.785 20 0 0 0 20 0 0 0
|
||||
3 1 1 2 4 5860.43 6.785 20 0 0 0 40 0 0 0
|
||||
4 1 1 3 5 5860.43 6.785 20 0 0 0 60 0 0 0
|
||||
5 1 1 4 6 5860.43 6.785 20 0 0 0 80 0 0 0
|
||||
6 1 1 5 7 5860.43 6.785 20 0 0 0 100 0 0 0
|
||||
7 1 1 6 8 5860.43 6.785 20 0 0 0 120 0 0 0
|
||||
8 1 1 7 9 5860.43 6.785 20 0 0 0 140 0 0 0
|
||||
9 1 1 8 10 5860.43 6.785 20 0 0 0 160 0 0 0
|
||||
10 1 1 9 11 5860.43 6.785 20 0 0 0 180 0 0 0
|
||||
11 1 1 10 1 5860.43 6.785 20 0 0 0 200 0 0 0
|
||||
12 2 1 22 13 5860.43 6.785 20 0 16.6992 0 0 0 0 0
|
||||
13 2 1 12 14 5860.43 6.785 20 0 16.6992 0 20 0 0 0
|
||||
14 2 1 13 15 5860.43 6.785 20 0 16.6992 0 40 0 0 0
|
||||
15 2 1 14 16 5860.43 6.785 20 0 16.6992 0 60 0 0 0
|
||||
16 2 1 15 17 5860.43 6.785 20 0 16.6992 0 80 0 0 0
|
||||
17 2 1 16 18 5860.43 6.785 20 0 16.6992 0 100 0 0 0
|
||||
18 2 1 17 19 5860.43 6.785 20 0 16.6992 0 120 0 0 0
|
||||
19 2 1 18 20 5860.43 6.785 20 0 16.6992 0 140 0 0 0
|
||||
20 2 1 19 21 5860.43 6.785 20 0 16.6992 0 160 0 0 0
|
||||
21 2 1 20 22 5860.43 6.785 20 0 16.6992 0 180 0 0 0
|
||||
22 2 1 21 12 5860.43 6.785 20 0 16.6992 0 200 0 0 0
|
||||
23 3 1 33 24 5860.43 6.785 20 0 8.3496 14.4619 0 0 0 0
|
||||
24 3 1 23 25 5860.43 6.785 20 0 8.3496 14.4619 20 0 0 0
|
||||
25 3 1 24 26 5860.43 6.785 20 0 8.3496 14.4619 40 0 0 0
|
||||
26 3 1 25 27 5860.43 6.785 20 0 8.3496 14.4619 60 0 0 0
|
||||
27 3 1 26 28 5860.43 6.785 20 0 8.3496 14.4619 80 0 0 0
|
||||
28 3 1 27 29 5860.43 6.785 20 0 8.3496 14.4619 100 0 0 0
|
||||
29 3 1 28 30 5860.43 6.785 20 0 8.3496 14.4619 120 0 0 0
|
||||
30 3 1 29 31 5860.43 6.785 20 0 8.3496 14.4619 140 0 0 0
|
||||
31 3 1 30 32 5860.43 6.785 20 0 8.3496 14.4619 160 0 0 0
|
||||
32 3 1 31 33 5860.43 6.785 20 0 8.3496 14.4619 180 0 0 0
|
||||
33 3 1 32 23 5860.43 6.785 20 0 8.3496 14.4619 200 0 0 0
|
||||
34 4 1 44 35 5860.43 6.785 20 0 -8.3496 14.4619 0 0 0 0
|
||||
35 4 1 34 36 5860.43 6.785 20 0 -8.3496 14.4619 20 0 0 0
|
||||
36 4 1 35 37 5860.43 6.785 20 0 -8.3496 14.4619 40 0 0 0
|
||||
37 4 1 36 38 5860.43 6.785 20 0 -8.3496 14.4619 60 0 0 0
|
||||
38 4 1 37 39 5860.43 6.785 20 0 -8.3496 14.4619 80 0 0 0
|
||||
39 4 1 38 40 5860.43 6.785 20 0 -8.3496 14.4619 100 0 0 0
|
||||
40 4 1 39 41 5860.43 6.785 20 0 -8.3496 14.4619 120 0 0 0
|
||||
41 4 1 40 42 5860.43 6.785 20 0 -8.3496 14.4619 140 0 0 0
|
||||
42 4 1 41 43 5860.43 6.785 20 0 -8.3496 14.4619 160 0 0 0
|
||||
43 4 1 42 44 5860.43 6.785 20 0 -8.3496 14.4619 180 0 0 0
|
||||
44 4 1 43 34 5860.43 6.785 20 0 -8.3496 14.4619 200 0 0 0
|
||||
45 5 1 55 46 5860.43 6.785 20 0 -16.6992 0 0 0 0 0
|
||||
46 5 1 45 47 5860.43 6.785 20 0 -16.6992 0 20 0 0 0
|
||||
47 5 1 46 48 5860.43 6.785 20 0 -16.6992 0 40 0 0 0
|
||||
48 5 1 47 49 5860.43 6.785 20 0 -16.6992 0 60 0 0 0
|
||||
49 5 1 48 50 5860.43 6.785 20 0 -16.6992 0 80 0 0 0
|
||||
50 5 1 49 51 5860.43 6.785 20 0 -16.6992 0 100 0 0 0
|
||||
51 5 1 50 52 5860.43 6.785 20 0 -16.6992 0 120 0 0 0
|
||||
52 5 1 51 53 5860.43 6.785 20 0 -16.6992 0 140 0 0 0
|
||||
53 5 1 52 54 5860.43 6.785 20 0 -16.6992 0 160 0 0 0
|
||||
54 5 1 53 55 5860.43 6.785 20 0 -16.6992 0 180 0 0 0
|
||||
55 5 1 54 45 5860.43 6.785 20 0 -16.6992 0 200 0 0 0
|
||||
56 6 1 66 57 5860.43 6.785 20 0 -8.3496 -14.4619 0 0 0 0
|
||||
57 6 1 56 58 5860.43 6.785 20 0 -8.3496 -14.4619 20 0 0 0
|
||||
58 6 1 57 59 5860.43 6.785 20 0 -8.3496 -14.4619 40 0 0 0
|
||||
59 6 1 58 60 5860.43 6.785 20 0 -8.3496 -14.4619 60 0 0 0
|
||||
60 6 1 59 61 5860.43 6.785 20 0 -8.3496 -14.4619 80 0 0 0
|
||||
61 6 1 60 62 5860.43 6.785 20 0 -8.3496 -14.4619 100 0 0 0
|
||||
62 6 1 61 63 5860.43 6.785 20 0 -8.3496 -14.4619 120 0 0 0
|
||||
63 6 1 62 64 5860.43 6.785 20 0 -8.3496 -14.4619 140 0 0 0
|
||||
64 6 1 63 65 5860.43 6.785 20 0 -8.3496 -14.4619 160 0 0 0
|
||||
65 6 1 64 66 5860.43 6.785 20 0 -8.3496 -14.4619 180 0 0 0
|
||||
66 6 1 65 56 5860.43 6.785 20 0 -8.3496 -14.4619 200 0 0 0
|
||||
67 7 1 77 68 5860.43 6.785 20 0 8.3496 -14.4619 0 0 0 0
|
||||
68 7 1 67 69 5860.43 6.785 20 0 8.3496 -14.4619 20 0 0 0
|
||||
69 7 1 68 70 5860.43 6.785 20 0 8.3496 -14.4619 40 0 0 0
|
||||
70 7 1 69 71 5860.43 6.785 20 0 8.3496 -14.4619 60 0 0 0
|
||||
71 7 1 70 72 5860.43 6.785 20 0 8.3496 -14.4619 80 0 0 0
|
||||
72 7 1 71 73 5860.43 6.785 20 0 8.3496 -14.4619 100 0 0 0
|
||||
73 7 1 72 74 5860.43 6.785 20 0 8.3496 -14.4619 120 0 0 0
|
||||
74 7 1 73 75 5860.43 6.785 20 0 8.3496 -14.4619 140 0 0 0
|
||||
75 7 1 74 76 5860.43 6.785 20 0 8.3496 -14.4619 160 0 0 0
|
||||
76 7 1 75 77 5860.43 6.785 20 0 8.3496 -14.4619 180 0 0 0
|
||||
77 7 1 76 67 5860.43 6.785 20 0 8.3496 -14.4619 200 0 0 0
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,29 +0,0 @@
|
||||
processors 1 1 *
|
||||
newton on
|
||||
units metal
|
||||
lattice sc 1.0
|
||||
boundary fs fs p
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 5 delay 0 check yes
|
||||
|
||||
atom_style mesont
|
||||
# cut, path, BendingMode, TPMType
|
||||
pair_style mesont/tpm 45.0 MESONT-TABTP_10_10.xrs 0 0
|
||||
read_data data.bundle
|
||||
pair_coeff * *
|
||||
|
||||
velocity all create 6000.0 2019
|
||||
timestep 0.005
|
||||
fix 1 all nve
|
||||
thermo 10
|
||||
reset_timestep 0
|
||||
|
||||
compute Es all mesont estretch
|
||||
compute Eb all mesont ebend
|
||||
compute Et all mesont etube
|
||||
compute B all property/atom buckling
|
||||
|
||||
thermo_style custom step time temp etotal ke pe c_Es c_Eb c_Et
|
||||
#dump out_dump all custom 50 dump.bundle id type x y z c_Es c_Eb c_Et c_B ix iy iz
|
||||
|
||||
run 100
|
||||
@ -1,28 +0,0 @@
|
||||
newton on
|
||||
units metal
|
||||
lattice sc 1.0
|
||||
boundary p p fs
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 5 delay 0 check yes
|
||||
|
||||
atom_style mesont
|
||||
# cut, path, BendingMode, TPMType
|
||||
pair_style mesont/tpm 30.0 MESONT-TABTP_10_10.xrs 1 0
|
||||
read_data data.film
|
||||
pair_coeff * *
|
||||
|
||||
velocity all create 600.0 2019
|
||||
timestep 0.01
|
||||
fix 1 all nve
|
||||
thermo 10
|
||||
reset_timestep 0
|
||||
|
||||
compute Es all mesont estretch
|
||||
compute Eb all mesont ebend
|
||||
compute Et all mesont etube
|
||||
compute B all property/atom buckling
|
||||
|
||||
thermo_style custom step time temp etotal ke pe c_Es c_Eb c_Et
|
||||
#dump out_dump all custom 10 dump.film id type x y z c_Es c_Eb c_Et c_B ix iy iz
|
||||
|
||||
run 20
|
||||
@ -1,90 +0,0 @@
|
||||
LAMMPS (5 May 2020)
|
||||
processors 1 1 *
|
||||
newton on
|
||||
units metal
|
||||
lattice sc 1.0
|
||||
Lattice spacing in x,y,z = 1 1 1
|
||||
boundary fs fs p
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 5 delay 0 check yes
|
||||
|
||||
atom_style mesont
|
||||
# cut, path, BendingMode, TPMType
|
||||
pair_style mesont/tpm 45.0 MESONT-TABTP_10_10.xrs 0 0
|
||||
read_data data.bundle
|
||||
orthogonal box = (-143.89 -143.89 0) to (143.89 143.89 220)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
77 atoms
|
||||
read_data CPU = 0.000613213 secs
|
||||
pair_coeff * *
|
||||
|
||||
velocity all create 6000.0 2019
|
||||
timestep 0.005
|
||||
fix 1 all nve
|
||||
thermo 10
|
||||
reset_timestep 0
|
||||
|
||||
compute Es all mesont estretch
|
||||
compute Eb all mesont ebend
|
||||
compute Et all mesont etube
|
||||
compute B all property/atom buckling
|
||||
|
||||
thermo_style custom step time temp etotal ke pe c_Es c_Eb c_Et
|
||||
#dump out_dump all custom 50 dump.bundle id type x y z c_Es c_Eb c_Et c_B ix iy iz
|
||||
|
||||
run 100
|
||||
Neighbor list info ...
|
||||
update every 5 steps, delay 0 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 46
|
||||
ghost atom cutoff = 46
|
||||
binsize = 23, bins = 7 7 10
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair mesont/tpm, perpetual
|
||||
attributes: full, newton on, ghost
|
||||
pair build: full/bin/ghost
|
||||
stencil: full/ghost/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 4.675 | 4.675 | 4.675 Mbytes
|
||||
Step Time Temp TotEng KinEng PotEng c_Es c_Eb c_Et
|
||||
0 0 6000 -201.86935 58.942626 -260.81198 0 0 -260.81198
|
||||
10 0.05 5114.1875 -201.86234 50.240607 -252.10295 4.8334861 2.3998206 -259.33626
|
||||
20 0.1 3437.2958 -201.8522 33.767207 -235.61941 11.42384 8.3426957 -255.38594
|
||||
30 0.15 2430.6571 -201.85242 23.878219 -225.73064 10.346152 14.72688 -250.80367
|
||||
40 0.2 2154.4755 -201.85683 21.165074 -223.0219 6.8146112 18.325709 -248.16222
|
||||
50 0.25 2021.7899 -201.85503 19.861601 -221.71663 9.2972022 17.644143 -248.65798
|
||||
60 0.3 2234.553 -201.85193 21.951737 -223.80367 13.541921 13.673721 -251.01931
|
||||
70 0.35 3099.6503 -201.85721 30.450255 -232.30747 11.833679 9.0583807 -253.19953
|
||||
80 0.4 3849.9855 -201.8635 37.821376 -239.68487 7.9899173 6.4332848 -254.10807
|
||||
90 0.45 3618.1311 -201.85967 35.543692 -237.40336 9.2616931 7.0452637 -253.71032
|
||||
100 0.5 2866.2722 -201.85273 28.157602 -230.01033 12.204916 10.284525 -252.49977
|
||||
Loop time of 0.419735 on 1 procs for 100 steps with 77 atoms
|
||||
|
||||
Performance: 102.922 ns/day, 0.233 hours/ns, 238.245 timesteps/s
|
||||
99.8% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.41922 | 0.41922 | 0.41922 | 0.0 | 99.88
|
||||
Neigh | 5.8174e-05 | 5.8174e-05 | 5.8174e-05 | 0.0 | 0.01
|
||||
Comm | 7.0333e-05 | 7.0333e-05 | 7.0333e-05 | 0.0 | 0.02
|
||||
Output | 0.00017667 | 0.00017667 | 0.00017667 | 0.0 | 0.04
|
||||
Modify | 0.00011945 | 0.00011945 | 0.00011945 | 0.0 | 0.03
|
||||
Other | | 8.702e-05 | | | 0.02
|
||||
|
||||
Nlocal: 77 ave 77 max 77 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 35 ave 35 max 35 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 2222 ave 2222 max 2222 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 2222
|
||||
Ave neighs/atom = 28.8571
|
||||
Neighbor list builds = 1
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:01
|
||||
@ -1,90 +0,0 @@
|
||||
LAMMPS (5 May 2020)
|
||||
processors 1 1 *
|
||||
newton on
|
||||
units metal
|
||||
lattice sc 1.0
|
||||
Lattice spacing in x,y,z = 1 1 1
|
||||
boundary fs fs p
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 5 delay 0 check yes
|
||||
|
||||
atom_style mesont
|
||||
# cut, path, BendingMode, TPMType
|
||||
pair_style mesont/tpm 45.0 MESONT-TABTP_10_10.xrs 0 0
|
||||
read_data data.bundle
|
||||
orthogonal box = (-143.89 -143.89 0) to (143.89 143.89 220)
|
||||
1 by 1 by 4 MPI processor grid
|
||||
reading atoms ...
|
||||
77 atoms
|
||||
read_data CPU = 0.000590563 secs
|
||||
pair_coeff * *
|
||||
|
||||
velocity all create 6000.0 2019
|
||||
timestep 0.005
|
||||
fix 1 all nve
|
||||
thermo 10
|
||||
reset_timestep 0
|
||||
|
||||
compute Es all mesont estretch
|
||||
compute Eb all mesont ebend
|
||||
compute Et all mesont etube
|
||||
compute B all property/atom buckling
|
||||
|
||||
thermo_style custom step time temp etotal ke pe c_Es c_Eb c_Et
|
||||
#dump out_dump all custom 50 dump.bundle id type x y z c_Es c_Eb c_Et c_B ix iy iz
|
||||
|
||||
run 100
|
||||
Neighbor list info ...
|
||||
update every 5 steps, delay 0 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 46
|
||||
ghost atom cutoff = 46
|
||||
binsize = 23, bins = 7 7 10
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair mesont/tpm, perpetual
|
||||
attributes: full, newton on, ghost
|
||||
pair build: full/bin/ghost
|
||||
stencil: full/ghost/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 4.671 | 4.671 | 4.671 Mbytes
|
||||
Step Time Temp TotEng KinEng PotEng c_Es c_Eb c_Et
|
||||
0 0 6000 -201.86935 58.942626 -260.81198 0 0 -260.81198
|
||||
10 0.05 5114.1875 -201.86234 50.240607 -252.10295 4.8334861 2.3998206 -259.33626
|
||||
20 0.1 3437.2958 -201.8522 33.767207 -235.61941 11.42384 8.3426957 -255.38594
|
||||
30 0.15 2430.6571 -201.85242 23.878219 -225.73064 10.346152 14.72688 -250.80367
|
||||
40 0.2 2154.4755 -201.85683 21.165074 -223.0219 6.8146112 18.325709 -248.16222
|
||||
50 0.25 2021.7899 -201.85503 19.861601 -221.71663 9.2972022 17.644143 -248.65798
|
||||
60 0.3 2234.553 -201.85193 21.951737 -223.80367 13.541921 13.673721 -251.01931
|
||||
70 0.35 3099.6503 -201.85721 30.450255 -232.30747 11.833679 9.0583807 -253.19953
|
||||
80 0.4 3849.9855 -201.8635 37.821376 -239.68487 7.9899173 6.4332848 -254.10807
|
||||
90 0.45 3618.1311 -201.85967 35.543692 -237.40336 9.2616931 7.0452637 -253.71032
|
||||
100 0.5 2866.2722 -201.85273 28.157602 -230.01033 12.204916 10.284525 -252.49977
|
||||
Loop time of 0.145545 on 4 procs for 100 steps with 77 atoms
|
||||
|
||||
Performance: 296.815 ns/day, 0.081 hours/ns, 687.071 timesteps/s
|
||||
95.9% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.046723 | 0.097529 | 0.14388 | 11.0 | 67.01
|
||||
Neigh | 2.5511e-05 | 2.8729e-05 | 3.171e-05 | 0.0 | 0.02
|
||||
Comm | 0.00058556 | 0.045174 | 0.098462 | 16.5 | 31.04
|
||||
Output | 0.0001483 | 0.0010182 | 0.002851 | 3.5 | 0.70
|
||||
Modify | 3.8147e-05 | 4.065e-05 | 4.4107e-05 | 0.0 | 0.03
|
||||
Other | | 0.001755 | | | 1.21
|
||||
|
||||
Nlocal: 19.25 ave 21 max 16 min
|
||||
Histogram: 1 0 0 0 0 0 1 0 0 2
|
||||
Nghost: 33.25 ave 40 max 28 min
|
||||
Histogram: 2 0 0 0 0 0 0 1 0 1
|
||||
Neighs: 0 ave 0 max 0 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 555.5 ave 606 max 460 min
|
||||
Histogram: 1 0 0 0 0 0 1 0 0 2
|
||||
|
||||
Total # of neighbors = 2222
|
||||
Ave neighs/atom = 28.8571
|
||||
Neighbor list builds = 1
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:01
|
||||
@ -1,81 +0,0 @@
|
||||
LAMMPS (5 May 2020)
|
||||
newton on
|
||||
units metal
|
||||
lattice sc 1.0
|
||||
Lattice spacing in x,y,z = 1 1 1
|
||||
boundary p p fs
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 5 delay 0 check yes
|
||||
|
||||
atom_style mesont
|
||||
# cut, path, BendingMode, TPMType
|
||||
pair_style mesont/tpm 30.0 MESONT-TABTP_10_10.xrs 1 0
|
||||
read_data data.film
|
||||
orthogonal box = (-2500 -2500 -300) to (2500 2500 402.42)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
79596 atoms
|
||||
read_data CPU = 0.0860629 secs
|
||||
pair_coeff * *
|
||||
|
||||
velocity all create 600.0 2019
|
||||
timestep 0.01
|
||||
fix 1 all nve
|
||||
thermo 10
|
||||
reset_timestep 0
|
||||
|
||||
compute Es all mesont estretch
|
||||
compute Eb all mesont ebend
|
||||
compute Et all mesont etube
|
||||
compute B all property/atom buckling
|
||||
|
||||
thermo_style custom step time temp etotal ke pe c_Es c_Eb c_Et
|
||||
#dump out_dump all custom 10 dump.film id type x y z c_Es c_Eb c_Et c_B ix iy iz
|
||||
|
||||
run 20
|
||||
Neighbor list info ...
|
||||
update every 5 steps, delay 0 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 31
|
||||
ghost atom cutoff = 31
|
||||
binsize = 15.5, bins = 323 323 26
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair mesont/tpm, perpetual
|
||||
attributes: full, newton on, ghost
|
||||
pair build: full/bin/ghost
|
||||
stencil: full/ghost/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 37.43 | 37.43 | 37.43 Mbytes
|
||||
Step Time Temp TotEng KinEng PotEng c_Es c_Eb c_Et
|
||||
0 0 600 1347.2158 6173.0767 -4825.8609 28.669574 21.29406 -4875.8245
|
||||
10 0.1 389.40755 1373.7864 4006.4045 -2632.6181 848.00267 1404.4323 -4885.053
|
||||
20 0.2 313.65714 1399.9427 3227.0494 -1827.1067 1201.1732 1882.1342 -4910.4141
|
||||
Loop time of 10.2438 on 1 procs for 20 steps with 79596 atoms
|
||||
|
||||
Performance: 1.687 ns/day, 14.228 hours/ns, 1.952 timesteps/s
|
||||
99.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 10.226 | 10.226 | 10.226 | 0.0 | 99.82
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.00081086 | 0.00081086 | 0.00081086 | 0.0 | 0.01
|
||||
Output | 0.00046563 | 0.00046563 | 0.00046563 | 0.0 | 0.00
|
||||
Modify | 0.015165 | 0.015165 | 0.015165 | 0.0 | 0.15
|
||||
Other | | 0.001869 | | | 0.02
|
||||
|
||||
Nlocal: 79596 ave 79596 max 79596 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 1879 ave 1879 max 1879 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 642270 ave 642270 max 642270 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 642270
|
||||
Ave neighs/atom = 8.06912
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:11
|
||||
@ -1,81 +0,0 @@
|
||||
LAMMPS (5 May 2020)
|
||||
newton on
|
||||
units metal
|
||||
lattice sc 1.0
|
||||
Lattice spacing in x,y,z = 1 1 1
|
||||
boundary p p fs
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 5 delay 0 check yes
|
||||
|
||||
atom_style mesont
|
||||
# cut, path, BendingMode, TPMType
|
||||
pair_style mesont/tpm 30.0 MESONT-TABTP_10_10.xrs 1 0
|
||||
read_data data.film
|
||||
orthogonal box = (-2500 -2500 -300) to (2500 2500 402.42)
|
||||
2 by 2 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
79596 atoms
|
||||
read_data CPU = 0.0704217 secs
|
||||
pair_coeff * *
|
||||
|
||||
velocity all create 600.0 2019
|
||||
timestep 0.01
|
||||
fix 1 all nve
|
||||
thermo 10
|
||||
reset_timestep 0
|
||||
|
||||
compute Es all mesont estretch
|
||||
compute Eb all mesont ebend
|
||||
compute Et all mesont etube
|
||||
compute B all property/atom buckling
|
||||
|
||||
thermo_style custom step time temp etotal ke pe c_Es c_Eb c_Et
|
||||
#dump out_dump all custom 10 dump.film id type x y z c_Es c_Eb c_Et c_B ix iy iz
|
||||
|
||||
run 20
|
||||
Neighbor list info ...
|
||||
update every 5 steps, delay 0 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 31
|
||||
ghost atom cutoff = 31
|
||||
binsize = 15.5, bins = 323 323 26
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair mesont/tpm, perpetual
|
||||
attributes: full, newton on, ghost
|
||||
pair build: full/bin/ghost
|
||||
stencil: full/ghost/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 12.81 | 12.81 | 12.83 Mbytes
|
||||
Step Time Temp TotEng KinEng PotEng c_Es c_Eb c_Et
|
||||
0 0 600 1347.2158 6173.0767 -4825.8609 28.669574 21.29406 -4875.8245
|
||||
10 0.1 389.40755 1373.7864 4006.4045 -2632.6181 848.00267 1404.4323 -4885.053
|
||||
20 0.2 313.65714 1399.9427 3227.0494 -1827.1067 1201.1732 1882.1342 -4910.4141
|
||||
Loop time of 3.67186 on 4 procs for 20 steps with 79596 atoms
|
||||
|
||||
Performance: 4.706 ns/day, 5.100 hours/ns, 5.447 timesteps/s
|
||||
95.8% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 2.7317 | 3.1286 | 3.6556 | 18.8 | 85.20
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.0036943 | 0.53094 | 0.92822 | 45.8 | 14.46
|
||||
Output | 0.00026512 | 0.00035298 | 0.00055647 | 0.0 | 0.01
|
||||
Modify | 0.010463 | 0.010884 | 0.011153 | 0.3 | 0.30
|
||||
Other | | 0.001109 | | | 0.03
|
||||
|
||||
Nlocal: 19899 ave 21951 max 18667 min
|
||||
Histogram: 1 1 0 1 0 0 0 0 0 1
|
||||
Nghost: 970.75 ave 1031 max 920 min
|
||||
Histogram: 1 0 0 1 1 0 0 0 0 1
|
||||
Neighs: 0 ave 0 max 0 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 160568 ave 181723 max 147382 min
|
||||
Histogram: 1 0 2 0 0 0 0 0 0 1
|
||||
|
||||
Total # of neighbors = 642270
|
||||
Ave neighs/atom = 8.06912
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:04
|
||||
24
src/atom.cpp
24
src/atom.cpp
@ -186,12 +186,6 @@ Atom::Atom(LAMMPS *_lmp) : Pointers(_lmp)
|
||||
cc = cc_flux = nullptr;
|
||||
edpd_temp = edpd_flux = edpd_cv = nullptr;
|
||||
|
||||
// MESONT package
|
||||
|
||||
length = nullptr;
|
||||
buckling = nullptr;
|
||||
bond_nt = nullptr;
|
||||
|
||||
// MACHDYN package
|
||||
|
||||
contact_radius = nullptr;
|
||||
@ -527,12 +521,6 @@ void Atom::peratom_create()
|
||||
add_peratom("cc",&cc,DOUBLE,1);
|
||||
add_peratom("cc_flux",&cc_flux,DOUBLE,1,1); // set per-thread flag
|
||||
|
||||
// MESONT package
|
||||
|
||||
add_peratom("length",&length,DOUBLE,0);
|
||||
add_peratom("buckling",&buckling,INT,0);
|
||||
add_peratom("bond_nt",&bond_nt,tagintsize,2);
|
||||
|
||||
// SPH package
|
||||
|
||||
add_peratom("rho",&rho,DOUBLE,0);
|
||||
@ -2919,12 +2907,6 @@ void *Atom::extract(const char *name)
|
||||
if (strcmp(name,"cv") == 0) return (void *) cv;
|
||||
if (strcmp(name,"vest") == 0) return (void *) vest;
|
||||
|
||||
// MESONT package
|
||||
|
||||
if (strcmp(name,"length") == 0) return (void *) length;
|
||||
if (strcmp(name,"buckling") == 0) return (void *) buckling;
|
||||
if (strcmp(name,"bond_nt") == 0) return (void *) bond_nt;
|
||||
|
||||
// MACHDYN package
|
||||
|
||||
if (strcmp(name, "contact_radius") == 0) return (void *) contact_radius;
|
||||
@ -3043,12 +3025,6 @@ int Atom::extract_datatype(const char *name)
|
||||
if (strcmp(name,"cv") == 0) return LAMMPS_DOUBLE;
|
||||
if (strcmp(name,"vest") == 0) return LAMMPS_DOUBLE_2D;
|
||||
|
||||
// MESONT package
|
||||
|
||||
if (strcmp(name,"length") == 0) return LAMMPS_DOUBLE;
|
||||
if (strcmp(name,"buckling") == 0) return LAMMPS_INT;
|
||||
if (strcmp(name,"bond_nt") == 0) return LAMMPS_TAGINT_2D;
|
||||
|
||||
// MACHDYN package
|
||||
|
||||
if (strcmp(name, "contact_radius") == 0) return LAMMPS_DOUBLE;
|
||||
|
||||
@ -145,12 +145,6 @@ class Atom : protected Pointers {
|
||||
double *edpd_cv; // heat capacity
|
||||
int cc_species;
|
||||
|
||||
// MESONT package
|
||||
|
||||
double *length;
|
||||
int *buckling;
|
||||
tagint **bond_nt;
|
||||
|
||||
// MACHDYN package
|
||||
|
||||
double *contact_radius;
|
||||
|
||||
@ -33,7 +33,6 @@ lmp2arc convert LAMMPS output to Accelrys Insight format
|
||||
lmp2cfg convert LAMMPS output to CFG files for AtomEye viz
|
||||
magic patterns to detect LAMMPS files with the file(1) command
|
||||
matlab MatLab scripts for post-processing LAMMPS output
|
||||
mesont Tools for use with the MESONT package
|
||||
micelle2d create a data file of small lipid chains in solvent
|
||||
moltemplate Instructions for installing the Moltemplate builder program
|
||||
msi2lmp use Accelrys Insight code to setup LAMMPS input
|
||||
|
||||
@ -1,58 +0,0 @@
|
||||
=== MESONT tools ===
|
||||
===============================
|
||||
|
||||
The programs in this folder can be used to analyze the
|
||||
output of simulations using the CNT mesoscopic force
|
||||
field (MESONT).
|
||||
|
||||
dump2vtk.cpp converts output written in *.dump format (the
|
||||
sequence of columns must be "ATOMS id type x y z Es Eb Et
|
||||
Ek ix iy iz", the same as in the examples at examples\USER\mesont)
|
||||
into VTK format that can be visualized as a set of tubes in
|
||||
Paraview (or other packages). The executable takes 3 parameters:
|
||||
system.init - an input file with information about connections
|
||||
between cnt nodes, config.dump - LAMMPS output with snapshots,
|
||||
out - output folder for writing VTK files (must exist).
|
||||
|
||||
Code TMDPotGen is designed to generate ASCII text files TPMSSTP.xrs
|
||||
and TPMA.xrs containing tabulated tubular potentials for
|
||||
single-walled CNTs with a given chirality (m,n). The input
|
||||
parameters for the code must be provided in the form of an ASCII
|
||||
text file TMDPotGen.xdt. The output of the code are files TPMSSTP.xrs
|
||||
and TPMA.xrs. All parameters in the tables are given in metal units.
|
||||
The generation of the tables takes approximately 4 hours.
|
||||
|
||||
Code TMDGen is designed to generate initial samples composed of
|
||||
straight and dispersed nanotubes of given chirality and length at
|
||||
a given material density. In the generated samples, nanotubes are
|
||||
distributed with random positions and orientations. Both periodic
|
||||
and free boundary conditions are available along each axis of the
|
||||
system. The input parameters for the code must be provided in form
|
||||
of an ASCII text file TMDGen.xdt and include the following:
|
||||
LS0: sample size along z- and y-directions (A)
|
||||
HS0: sample size along z-direction (A)
|
||||
DS0: material density (g/cm^3)
|
||||
BC_X0: Type of boundary conditions along x-direction (0, Free; 1, Periodic)
|
||||
BC_Y0: Type of boundary conditions along y-direction (0, Free; 1, Periodic)
|
||||
BC_Z0: Type of boundary conditions along z-direction (0, Free; 1, Periodic)
|
||||
ChiIndM: First chirality index of nanotubes
|
||||
ChiIndN: Second chirality index of nanotubes
|
||||
LT0: Nanotube length (A)
|
||||
SegType: Parameter that defines how a nanotubes will be divided into
|
||||
segments(0, NSeg0 will be used; 1, LSeg0 will be used)
|
||||
NSeg0: Number of segments in every nanotube. Used if SegType = 0. Then
|
||||
LSeg0 = LT0 / NSeg0
|
||||
LSeg0: Length of segments in every nanotube. Used if SegType = 1. Then
|
||||
NSeg0 = [ LT0 / LSeg0 ]
|
||||
DeltaT: Minimum gap between nanotube walls in the generated sample (A)
|
||||
NAmax: Maximum number of attempts to add new nanotube to the sample
|
||||
GeomPrec: Precision of calculations (dimensionless).
|
||||
The output of the code is an ASCII text file TMDSample.init written in the
|
||||
LAMMPS format compatible with cnt atomic style. All parameters in the sample
|
||||
files generated with TMDGen are given in metal units.
|
||||
|
||||
|
||||
This packages were created by Maxim Shugaev (mvs9t@virginia.edu)
|
||||
at the University of Virginia and by Alexey N. Volkov (avolkov1@ua.edu)
|
||||
at the University of Alabama.
|
||||
|
||||
@ -1,33 +0,0 @@
|
||||
#---------------------------------------------------------------------------------------------------
|
||||
#
|
||||
# This is Makefile for builing the executable TMDGen
|
||||
#
|
||||
# Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, 2020, Version 13.00
|
||||
#
|
||||
#---------------------------------------------------------------------------------------------------
|
||||
EXEPATH = .
|
||||
|
||||
F90 = ifort
|
||||
F90FLAGS = -O3 -ipo
|
||||
LDFLAGS =
|
||||
|
||||
OBJS = TPMLib.o TPMGeom.o TMDGenData.o TMDGen3D.o TMDGen.o
|
||||
EXE = $(EXEPATH)/TMDGen
|
||||
|
||||
# compile and load
|
||||
default:
|
||||
@echo " "
|
||||
@echo "Compiling Code of Program TMDGen"
|
||||
@echo "FORTRAN 90"
|
||||
$(MAKE) $(EXE)
|
||||
|
||||
$(EXE): $(OBJS)
|
||||
$(F90) $(F90FLAGS) $(LDFLAGS) -o $(EXE) $(OBJS)
|
||||
|
||||
.SUFFIXES: .f90 .o
|
||||
.f90.o:
|
||||
$(F90) $(F90FLAGS) -c $*.f90
|
||||
|
||||
clean:
|
||||
rm -f *.o
|
||||
|
||||
@ -1,267 +0,0 @@
|
||||
program TMDGen !************************************************************************************
|
||||
!
|
||||
! Stand-alone generator of 3D CNT samples.
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, 2020, Version 13.00
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
use TMDGen3D
|
||||
|
||||
implicit none
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Global variables
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4 :: Nseg, Nnode
|
||||
real*8 :: DS00
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Body
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
print *, 'TMD generator of 3D CNT samples, v. 13.00'
|
||||
print '(a34,a,i10)', 'Maximum number of nanotubes', ' : ', MAX_TUBE
|
||||
|
||||
call SetRandomSeed ()
|
||||
|
||||
! Reading and printing of governing parameters
|
||||
call LoadGoverningParameters ()
|
||||
call PrintGoverningParameters ()
|
||||
|
||||
! Here we calculate the radius of nanotubes
|
||||
RT0 = TPBA * sqrt ( 3.0d+00 * ( ChiIndM * ChiIndM + ChiIndN * ChiIndN + ChiIndM * ChiIndN ) ) / M_2PI;
|
||||
|
||||
! Here we calculate parameters of the desired sample
|
||||
call InitSample ()
|
||||
DS0 = DS0 * ( K_MDDU / 1.0d+03 )
|
||||
call PrintSampleParameters ( 'Desired' )
|
||||
DS00 = DS0
|
||||
DS0 = DS0 / ( K_MDDU / 1.0d+03 )
|
||||
|
||||
call Generator3D ()
|
||||
|
||||
! Here we write the major output file with the sample
|
||||
!call WriteOutputFile_old_format ()
|
||||
!call WriteOutputFile ()
|
||||
|
||||
! Here we write an auxiliary Tecplot file to visualize the initial sample
|
||||
!PrintTecplotFile ()
|
||||
call WriteLAMMPSFile()
|
||||
|
||||
! Here we print parameters of the final sample
|
||||
call PrintSampleParameters ( 'Final' )
|
||||
print '(a34,a,f15.4,a)', 'Nanotube radius ', ' : ', RT0, ' a'
|
||||
print '(a34,a,f15.4,a)', 'Nanotube length ', ' : ', LT0, ' a'
|
||||
print '(a34,a,f15.4,a)', 'Nanotube mass ', ' : ', M_2PI * RT0 * LT0 * TPBM * TPBD, ' Da'
|
||||
if ( SegType == 0 ) then
|
||||
LSeg0 = LT0 / NSeg0
|
||||
else
|
||||
NSeg0 = int ( LT0 / LSeg0 ) + 1
|
||||
LSeg0 = LT0 / NSeg0
|
||||
end if
|
||||
print '(a34,a,f15.4,a)', 'Nanotube segment length ', ' : ', LSeg0, ' a'
|
||||
print '(a34,a,f15.4,a)', 'Nanotube segment mass ', ' : ', M_2PI * RT0 * LSeg0 * TPBM * TPBD, ' Da'
|
||||
print '(a34,a,f15.4)', 'Desired / Real densities ', ' : ', DS00 / DS0
|
||||
print '(a34,a,i10)', 'Real number of tubes', ' : ', NT
|
||||
print '(a34,a,i10)', 'Real number of segments', ' : ', Nseg
|
||||
print '(a34,a,i10)', 'Real number of nodes', ' : ', Nnode
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
subroutine DiscretizeTube ( X0, DL, NS, i ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function calculaats geometrical parameters that are necessary to represent straight
|
||||
! tube i as a sequence of segments.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2), intent(out) :: X0
|
||||
real*8, intent(out) :: DL
|
||||
integer*4, intent(out) :: NS
|
||||
integer*4, intent(in) :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2) :: X1
|
||||
!-------------------------------------------------------------------------------------------
|
||||
call GetTubeEnds ( X0, X1, i )
|
||||
if ( SegType == 0 ) then
|
||||
NS = NSeg0
|
||||
else
|
||||
NS = int ( LT(i) / LSeg0 ) + 1
|
||||
end if
|
||||
DL = LT(i) / NS
|
||||
end subroutine DiscretizeTube !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine WriteOutputFile_old_format () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function writes a dat file (version 2) with the initial nanotube sample.
|
||||
! This file is used by TMD/TMDMPI to start a new simulation.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: Fuid, i, j, NTS, Prop
|
||||
real*8 :: DL, L, L00, M00, I00, J00, C00, LL00, MM00, II00, JJ00, CC00
|
||||
real*8, dimension(0:2) :: X, X0
|
||||
logical*4 :: PrintNode
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Fuid = OpenFile ( 'TMDGen_old.dat', "wt", "" )
|
||||
write ( unit = Fuid, fmt = '(i12)' ) 3
|
||||
write ( unit = Fuid, fmt = '(2i4,4e20.12)' ) ChiIndM, ChiIndN, RT0, TPBA, TPBD, TPBM
|
||||
write ( unit = Fuid, fmt = '(3e20.12)' ) DomXmin, DomYmin, DomZmin
|
||||
write ( unit = Fuid, fmt = '(3e20.12)' ) DomXmax, DomYmax, DomZmax
|
||||
write ( unit = Fuid, fmt = '(3i12)' ) BC_X, BC_Y, BC_Z
|
||||
write ( unit = Fuid, fmt = '(i12)' ) NT
|
||||
Nseg = 0
|
||||
Nnode = 0
|
||||
do i = 0, NT - 1
|
||||
call DiscretizeTube ( X0, DL, NTS, i )
|
||||
L00 = LT(i) / NTS
|
||||
M00 = TubeMass ( i ) / NTS
|
||||
I00 = 0.0d+00
|
||||
J00 = M00 * sqr ( RT(i) )
|
||||
C00 = M00 * TubeSpecificHeat ( i )
|
||||
Nseg = Nseg + NTS
|
||||
write ( unit = Fuid, fmt = '(i12)' ) NTS + 1
|
||||
Nnode = Nnode + NTS + 1
|
||||
L = 0.0d+00
|
||||
do j = 0, NTS
|
||||
X = X0 + L * DT(i,0:2)
|
||||
MM00 = M00
|
||||
II00 = I00
|
||||
JJ00 = J00
|
||||
CC00 = C00
|
||||
LL00 = L00
|
||||
if ( j == 0 .or. j == NTS ) then
|
||||
MM00 = 0.5d+00 * M00
|
||||
II00 = 0.5d+00 * I00
|
||||
JJ00 = 0.5d+00 * J00
|
||||
CC00 = 0.5d+00 * C00
|
||||
end if
|
||||
if ( j == NTS ) LL00 = 0.0d+00
|
||||
Prop = 0
|
||||
write ( unit = Fuid, fmt = '(i2,6e20.12)' ) Prop, RT(0), LL00, MM00, II00, JJ00, CC00
|
||||
write ( unit = Fuid, fmt = '(6e20.12)' ) X, RT(i), 0.0d+00, 300.0d+00
|
||||
L = L + DL
|
||||
end do
|
||||
end do
|
||||
write ( unit = Fuid, fmt = '(i12)' ) 0
|
||||
write ( unit = Fuid, fmt = '(i12)' ) 0
|
||||
call CloseFile ( Fuid )
|
||||
end subroutine WriteOutputFile_old_format !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine WriteOutputFile () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function writes a dat file (version 2) with the initial nanotube sample.
|
||||
! This file is used by TMD/TMDMPI to start a new simulation.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: Fuid, i, j, NTS
|
||||
real*8 :: DL, L, L00, M00, LL00, MM00
|
||||
real*8, dimension(0:2) :: X, X0
|
||||
logical*4 :: PrintNode
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Fuid = OpenFile ( 'TMDGen.dat', "wt", "" )
|
||||
write ( unit = Fuid, fmt = '(2i4,4e20.12)' ) ChiIndM, ChiIndN, RT0, TPBA, TPBD, TPBM
|
||||
write ( unit = Fuid, fmt = '(3e20.12)' ) DomXmin, DomYmin, DomZmin
|
||||
write ( unit = Fuid, fmt = '(3e20.12)' ) DomXmax, DomYmax, DomZmax
|
||||
write ( unit = Fuid, fmt = '(3i12)' ) BC_X, BC_Y, BC_Z
|
||||
write ( unit = Fuid, fmt = '(i12)' ) NT
|
||||
Nseg = 0
|
||||
Nnode = 0
|
||||
do i = 0, NT - 1
|
||||
call DiscretizeTube ( X0, DL, NTS, i )
|
||||
L00 = LT(i) / NTS
|
||||
M00 = TubeMass ( i ) / NTS
|
||||
Nseg = Nseg + NTS
|
||||
write ( unit = Fuid, fmt = '(i12)' ) NTS + 1
|
||||
Nnode = Nnode + NTS + 1
|
||||
L = 0.0d+00
|
||||
do j = 0, NTS
|
||||
X = X0 + L * DT(i,0:2)
|
||||
MM00 = M00
|
||||
LL00 = L00
|
||||
if ( j == 0 .or. j == NTS ) MM00 = 0.5d+00 * M00
|
||||
if ( j == NTS ) LL00 = 0.0d+00
|
||||
write ( unit = Fuid, fmt = '(5e20.12)' ) X, LL00, MM00
|
||||
L = L + DL
|
||||
end do
|
||||
end do
|
||||
call CloseFile ( Fuid )
|
||||
end subroutine WriteOutputFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine PrintTecplotFile () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function prints Tecplot file to visualize the generated sample
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: Fuid, i
|
||||
real*8 :: LT2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Fuid = OpenFile ( 'TMDGen.plt', "wt", "" )
|
||||
write ( unit = Fuid, fmt = '(a)' ) 'VARIABLES="X" "Y" "Z"'
|
||||
do i = 0, NT - 1
|
||||
write ( unit = Fuid, fmt = '(a,i,a)' ) 'ZONE T="T', i, '"'
|
||||
LT2 = 0.5d+00 * LT(i)
|
||||
write ( unit = Fuid, fmt = '(3e20.12)' ) CT(i,0:2) - LT2 * DT(i,0:2)
|
||||
write ( unit = Fuid, fmt = '(3e20.12)' ) CT(i,0:2) + LT2 * DT(i,0:2)
|
||||
end do
|
||||
call CloseFile ( Fuid )
|
||||
end subroutine PrintTecplotFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
|
||||
subroutine WriteLAMMPSFile () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function writes a dat file (version 2) with the initial nanotube sample.
|
||||
! This file is used by TMD/TMDMPI to start a new simulation.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: file_id, i, j, NTS, node_id, b1, b2
|
||||
real*8 :: DL, L, L00, M00, LL00, MM00
|
||||
real*8, dimension(0:2) :: X, X0
|
||||
logical*4 :: PrintNode
|
||||
!-------------------------------------------------------------------------------------------
|
||||
open(newunit = file_id, file = 'TMDSample.init')
|
||||
write(file_id,*)
|
||||
write(file_id,*)
|
||||
!count the number of nodes and segments
|
||||
Nseg = 0
|
||||
Nnode = 0
|
||||
do i = 0, NT - 1
|
||||
call DiscretizeTube (X0, DL, NTS, i)
|
||||
Nseg = Nseg + NTS
|
||||
Nnode = Nnode + NTS + 1
|
||||
enddo
|
||||
write(file_id,'(i9,a)') Nnode, " atoms"
|
||||
write(file_id,*)
|
||||
write(file_id,*) "1 atom types"
|
||||
write(file_id,*)
|
||||
write(file_id,'(2e20.12,2a)') DomXmin, DomXmax, " xlo xhi"
|
||||
write(file_id,'(2e20.12,2a)') DomYmin, DomYmax, " ylo yhi"
|
||||
write(file_id,'(2e20.12,2a)') DomZmin, DomZmax, " zlo zhi"
|
||||
write(file_id,*)
|
||||
write(file_id,*) "Masses"
|
||||
write(file_id,*)
|
||||
write(file_id,*) "1 1.0"
|
||||
write(file_id,*)
|
||||
write(file_id,*) "Atoms"
|
||||
write(file_id,*)
|
||||
|
||||
node_id = 1
|
||||
do i = 0, NT - 1
|
||||
call DiscretizeTube(X0, DL, NTS, i)
|
||||
L00 = LT(i) / NTS
|
||||
M00 = TubeMass (i) / NTS
|
||||
|
||||
b1 = -1
|
||||
L = 0.0d+00
|
||||
do j = 0, NTS
|
||||
b2 = node_id + 1
|
||||
if (j == NTS) b2 = -1
|
||||
MM00 = M00
|
||||
LL00 = L00
|
||||
if (j == 0 .or. j == NTS) MM00 = 0.5d+00 * M00
|
||||
if (j == NTS) LL00 = 0.0d+00
|
||||
X = X0 + L * DT(i,0:2)
|
||||
write(file_id,'(2i9,a,2i9,3e14.7,a,3e20.12,a)') node_id, i, " 1 ", b1, b2, MM00, RT(i), LL00, " 0 ", X, " 0 0 0"
|
||||
b1 = node_id
|
||||
node_id = node_id + 1
|
||||
L = L + DL
|
||||
enddo
|
||||
enddo
|
||||
close(file_id)
|
||||
end subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end program TMDGen !********************************************************************************
|
||||
@ -1,15 +0,0 @@
|
||||
0.400000000000E+04 : LS0, A
|
||||
0.400000000000E+04 : HS0, A
|
||||
0.010000000000E+00 : DS0, Density g/cm^3
|
||||
1 : BC_X0, periodic along X
|
||||
1 : BC_Y0, periodic along Y
|
||||
0 : BC_Z0, periodic along Z
|
||||
10 : ChiIndM, tube chirality M
|
||||
10 : ChiIndN, tube chirality N
|
||||
0.200000000000E+04 : LT0, A
|
||||
0 : SegType
|
||||
100 : NSeg0
|
||||
0.200000000000E+02 : LSeg0
|
||||
0.500000000000E+01 : DeltaT, A
|
||||
1000000 : NAmax
|
||||
0.100000000000E-06 : GeomPrec
|
||||
@ -1,231 +0,0 @@
|
||||
module TMDGen3D !***********************************************************************************
|
||||
!
|
||||
! Generator of 3D CNT samples for TPM force field
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 13.00, 2020
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
use TMDGenData
|
||||
|
||||
implicit none
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
real*8 function MinimalDistance3D ( S1, S2, H, cosA, P, Q ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function returns the minimum distance between two line segments in 3D
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(out) :: S1, S2
|
||||
real*8, intent(in) :: H, cosA
|
||||
real*8, dimension(0:1), intent(in) :: P, Q
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: H2, cosA2, D
|
||||
real*8, dimension(0:1) :: P1, Q1
|
||||
integer*4, dimension(0:1,0:1) :: KA
|
||||
integer*4 :: i, j, K
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( ( P(0) * P(1) .le. 0.0d+00 ) .and. ( Q(0) * Q(1) .le. 0.0d+00 ) ) then
|
||||
MinimalDistance3D = H
|
||||
S1 = 0.5d+00 * ( P(0) + P(1) )
|
||||
S2 = 0.5d+00 * ( Q(0) + Q(1) )
|
||||
return
|
||||
end if
|
||||
do i = 0, 1
|
||||
P1(i) = P(i) * cosA
|
||||
Q1(i) = Q(i) * cosA
|
||||
end do
|
||||
KA = 1
|
||||
K = 0
|
||||
do i = 0, 1
|
||||
if ( ( Q1(i) .ge. P(0) ) .and. ( Q1(i) .le. P(1) ) ) then
|
||||
D = sqr ( Q(i) )
|
||||
if ( K == 0 ) then
|
||||
MinimalDistance3D = D
|
||||
S1 = Q1(i)
|
||||
S2 = Q(i)
|
||||
K = 1
|
||||
else if ( D < MinimalDistance3D ) then
|
||||
MinimalDistance3D = D
|
||||
S1 = Q1(i)
|
||||
S2 = Q(i)
|
||||
end if
|
||||
KA(0,i) = 0
|
||||
KA(1,i) = 0
|
||||
end if
|
||||
if ( ( P1(i) .ge. Q(0) ) .and. ( P1(i) .le. Q(1) ) ) then
|
||||
D = sqr ( P(i) )
|
||||
if ( K == 0 ) then
|
||||
MinimalDistance3D = D
|
||||
S1 = P(i)
|
||||
S2 = P1(i)
|
||||
K = 1
|
||||
else if ( D < MinimalDistance3D ) then
|
||||
MinimalDistance3D = D
|
||||
S1 = P(i)
|
||||
S2 = P1(i)
|
||||
end if
|
||||
KA(i,0) = 0
|
||||
KA(i,1) = 0
|
||||
end if
|
||||
end do
|
||||
H2 = sqr ( H )
|
||||
cosA2 = 2.0d+00 * cosA
|
||||
if ( K == 1 ) MinimalDistance3D = H2 + MinimalDistance3D * ( 1.0d+00 - sqr ( cosA ) )
|
||||
do i = 0, 1
|
||||
do j = 0, 1
|
||||
if ( KA(i,j) == 1 ) then
|
||||
D = H2 + sqr ( P(i) ) + sqr ( Q(j) ) - P(i) * Q(j) * cosA2
|
||||
if ( K == 0 ) then
|
||||
MinimalDistance3D = D
|
||||
S1 = P(i)
|
||||
S2 = Q(j)
|
||||
K = 1
|
||||
else if ( D < MinimalDistance3D ) then
|
||||
MinimalDistance3D = D
|
||||
S1 = P(i)
|
||||
S2 = Q(j)
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
MinimalDistance3D = dsqrt ( MinimalDistance3D )
|
||||
end function MinimalDistance3D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine RandTube3D ( X, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This subroutine generates a random tube in an isotropic 3D sample
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2), intent(out) :: X, L
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: CT, ST, E
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( BC_X0 == 0 ) then
|
||||
X(0)= LS0 * randnumber ()
|
||||
else
|
||||
X(0)= LS0 * ( 0.5d+00 - 1.0d+00 * randnumber () )
|
||||
end if
|
||||
if ( BC_Y0 == 0 ) then
|
||||
X(1)= LS0 * randnumber ()
|
||||
else
|
||||
X(1)= LS0 * ( 0.5d+00 - 1.0d+00 * randnumber () )
|
||||
end if
|
||||
if ( BC_Z0 == 0 ) then
|
||||
X(2)= HS0 *randnumber ()
|
||||
else
|
||||
X(2)= HS0 * ( 0.5d+00 - 1.0d+00 * randnumber () )
|
||||
end if
|
||||
CT = 1.0d+00 - 2.0d+00 * randnumber ()
|
||||
ST = sqrt ( 1.0d+00 - sqr ( CT ) )
|
||||
E = M_2PI * randnumber ()
|
||||
L(0)= CT
|
||||
L(1)= ST * cos ( E )
|
||||
L(2)= ST * sin ( E )
|
||||
end subroutine RandTube3D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
logical function AddTubeToSample3D ( MS ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function adds the last generated tube to the existing sample, if possible.
|
||||
! In a case of periodic boundaries, this version is valid only f the tube length is smaller
|
||||
! than the half of the sample.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(inout) :: MS
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: i, m
|
||||
real*8 :: Dmin, LT2, H, cosA, D1, D2, S1, S2
|
||||
real*8, dimension(0:2) :: X, L12
|
||||
real*8, dimension(0:1) :: P, Q
|
||||
!-------------------------------------------------------------------------------------------
|
||||
AddTubeToSample3D = .false.
|
||||
if ( .not. IsTubeInside ( NT ) ) return
|
||||
|
||||
LT2 = 0.5d+00 * LT(NT)
|
||||
do m = 0, NT - 1
|
||||
X = CT(NT,0:2)
|
||||
if ( LineLine ( H, cosA, D1, D2, L12, X, DT(NT,0:2), CT(m,0:2), DT(m,0:2), GeomPrec ) == MD_LINES_NONPAR ) then
|
||||
P(0) = D1 - LT2
|
||||
P(1) = D1 + LT2
|
||||
Q(0) = D2 - 0.5d+00 * LT(m)
|
||||
Q(1) = D2 + 0.5d+00 * LT(m)
|
||||
Dmin = MinimalDistance3D ( S1, S2, H, cosA, P, Q )
|
||||
else
|
||||
call LinePoint ( H, L12, CT(m,0:2), DT(m,0:2), X )
|
||||
L12 = L12 - X
|
||||
call ApplyPeriodicBC ( L12 )
|
||||
Dmin = S_V3norm3 ( L12 )
|
||||
end if
|
||||
if ( Dmin < RT(NT) + RT(m) + DeltaT ) return
|
||||
end do
|
||||
|
||||
MS = MS + TubeMass ( NT )
|
||||
NT = NT + 1
|
||||
AddTubeToSample3D = .true.
|
||||
end function AddTubeToSample3D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine Generator3D () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This subroutine implements the whole fgenerator of 3D samples
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: NA, NT0
|
||||
real*8 :: MS
|
||||
real*8 :: X1, X2, Y1, Y2, Z1, Z2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
NT = 0
|
||||
MS = 0.0d+00
|
||||
NT0 = int ( MS0 / ( M_2PI * RT0 * LT0 * TPBM * TPBD ) )
|
||||
do
|
||||
if ( NT == MAX_TUBE ) then
|
||||
print *, 'Error in [Generator3D]: MAX_TUBE is too small'
|
||||
stop
|
||||
end if
|
||||
if ( MS .ge. MS0 ) exit
|
||||
NA = 0
|
||||
! Trying to add the tube to the sample
|
||||
! The maximal number of attempts is equal to NAmax
|
||||
RT(NT) = RT0
|
||||
LT(NT) = LT0
|
||||
do
|
||||
if ( NA == NAmax ) exit
|
||||
call RandTube3D ( CT(NT,0:2), DT(NT,0:2) )
|
||||
if ( AddTubeToSample3D ( MS ) ) then
|
||||
print '(a,i10,a,i10,a,i10)', 'Tube ', NT, '(Appr.', NT0, ' total): Attempt ', NA
|
||||
if ( BC_X0 == 0 ) then
|
||||
X1 = CT(NT,0) - 0.5d+00 * LT(NT) * DT(NT,0)
|
||||
X2 = CT(NT,0) + 0.5d+00 * LT(NT) * DT(NT,0)
|
||||
if ( DomXmin > X1 ) DomXmin = X1
|
||||
if ( DomXmin > X2 ) DomXmin = X2
|
||||
if ( DomXmax < X1 ) DomXmax = X1
|
||||
if ( DomXmax < X2 ) DomXmax = X2
|
||||
end if
|
||||
if ( BC_Y0 == 0 ) then
|
||||
Y1 = CT(NT,1) - 0.5d+00 * LT(NT) * DT(NT,1)
|
||||
Y2 = CT(NT,1) + 0.5d+00 * LT(NT) * DT(NT,1)
|
||||
if ( DomYmin > Y1 ) DomYmin = Y1
|
||||
if ( DomYmin > Y2 ) DomYmin = Y2
|
||||
if ( DomYmax < Y1 ) DomYmax = Y1
|
||||
if ( DomYmax < Y2 ) DomYmax = Y2
|
||||
end if
|
||||
if ( BC_Z0 == 0 ) then
|
||||
Z1 = CT(NT,2) - 0.5d+00 * LT(NT) * DT(NT,2)
|
||||
Z2 = CT(NT,2) + 0.5d+00 * LT(NT) * DT(NT,2)
|
||||
if ( DomZmin > Z1 ) DomZmin = Z1
|
||||
if ( DomZmin > Z2 ) DomZmin = Z2
|
||||
if ( DomZmax < Z1 ) DomZmax = Z1
|
||||
if ( DomZmax < Z2 ) DomZmax = Z2
|
||||
end if
|
||||
exit
|
||||
end if
|
||||
NA = NA + 1
|
||||
end do
|
||||
end do
|
||||
MS0 = MS
|
||||
|
||||
if ( BC_X0 == 0 ) DomLX = DomXmax - DomXmin
|
||||
if ( BC_Y0 == 0 ) DomLY = DomYmax - DomYmin
|
||||
if ( BC_Z0 == 0 ) DomLZ = DomZmax - DomZmin
|
||||
|
||||
VS0 = ( DomXmax - DomXmin ) * ( DomYmax - DomYmin ) * ( DomZmax - DomZmin )
|
||||
DS0 = MS0 / VS0 * ( K_MDDU / 1.0d+03 )
|
||||
end subroutine Generator3D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module TMDGen3D !*******************************************************************************
|
||||
@ -1,289 +0,0 @@
|
||||
module TMDGenData !*********************************************************************************
|
||||
!
|
||||
! Common data for TMDGen
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 13.00, 2020
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
use TPMGeom
|
||||
|
||||
implicit none
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Constants
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4, parameter :: MAX_TUBE = 1000000 ! Maximum number of tubes in 3D
|
||||
|
||||
real*8, parameter :: K_MDDU = K_MDMU / K_MDLU / K_MDLU / K_MDLU ! MD density unit (kg/m**3)
|
||||
|
||||
!
|
||||
! These parameters are specific for carbon nanotubes and taken from module TubePotBase
|
||||
!
|
||||
|
||||
real*8, parameter :: TPbConstD = 5.196152422706632d+00 ! = 3.0**1.5
|
||||
! Mass of C atom
|
||||
real*8, parameter :: TPBM = 12.0107d+00 ! (a.m.u.)
|
||||
! Lattice parameter and numerical density of atoms for a graphene sheet, see Dresselhaus et al, Carbon 33(7), 1995
|
||||
real*8, parameter :: TPBA = 1.421d+00 ! (Angstrom)
|
||||
real*8, parameter :: TPBD = 4.0d+00 / ( TPBConstD * TPBA * TPBA ) ! (1/Angstrom^2)
|
||||
! Specific heat of carbon nanotubes
|
||||
real*8, parameter :: TPBSH = 600.0d+00 / K_MDCU ! (eV/(Da*K))
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Governing parameters
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
! Parameters of the sample
|
||||
|
||||
real*8 :: LS0 = 4000.0 ! Sample size in x, y-directions (Angstrom)
|
||||
real*8 :: HS0 = 4000.0 ! Sample size in z-direction (Angstrom)
|
||||
real*8 :: DS0 = 0.01 ! Density (g/cm**3)
|
||||
integer*4 :: BC_X0 = 1 ! Boundary conditions in x-direction: 0, free; 1, periodic
|
||||
integer*4 :: BC_Y0 = 1 ! Boundary conditions in y-direction: 0, free; 1, periodic
|
||||
integer*4 :: BC_Z0 = 1 ! Boundary conditions in z-direction: 0, free; 1, periodic
|
||||
|
||||
! Parameters of tubes
|
||||
integer*4 :: ChiIndM = 10 ! Chirality index m of nanotubes
|
||||
integer*4 :: ChiIndN = 10 ! Chirality index n of nanotubes
|
||||
real*8 :: LT0 = 2000.0 ! Characterstic length of tubes (Angstrom)
|
||||
integer*4 :: SegType = 0 ! 0, number of segments per tube is fixed
|
||||
! 1, rounded length of segments is fixed
|
||||
integer*4 :: NSeg0 = 100 ! Number of segments per tube
|
||||
real*8 :: LSeg0 = 20.0d+00 ! Length of the segment (Angstrom)
|
||||
|
||||
|
||||
! Parameters controlling the sample structure
|
||||
|
||||
real*8 :: DeltaT = 3.0 ! Minimal distance between tubes (Angstrom)
|
||||
integer*4 :: NAmax = 50000 ! Maximal number of attempts (for SampleType = 4 it is used as an input paramtere for number of tubes)
|
||||
real*8 :: GeomPrec = 1.0d-06 ! Geometrical precision
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Computed data
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8 :: RT0 = 6.785 ! Radius of tubes (Angstrom)
|
||||
|
||||
real*8 :: VS0 ! Desired volume of the sample, Angstrom**3
|
||||
real*8 :: MS0 ! Desired mass of the sample, Da (For SampleType = 4 it is the defined fixed mass- definition is given in TMDGen7T)
|
||||
|
||||
real*8 :: CTCD ! Center to center distance between any surrounding tube and center tube (used for SampleType == 4 only)
|
||||
|
||||
integer*4 :: NT ! Real number of tubes
|
||||
real*8, dimension(0:MAX_TUBE-1) :: RT ! Radii of tubes, Angstrom
|
||||
real*8, dimension(0:MAX_TUBE-1) :: LT ! Lengths of tubes, Angstrom
|
||||
real*8, dimension(0:MAX_TUBE-1,0:2) :: CT ! Coordinates of tubes' centers, Angstrom
|
||||
real*8, dimension(0:MAX_TUBE-1,0:2) :: DT ! Directions of tubes
|
||||
integer*4, dimension(0:MAX_TUBE-1) :: AT ! Parent axes of tubes. It is used only in GeneratorBundle ()
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Pseudo-random number generator
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8 function randnumber () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function returns a pseudo-random number with uniform distribution in [0,1]
|
||||
!-------------------------------------------------------------------------------------------
|
||||
call random_number ( randnumber )
|
||||
end function randnumber !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine SetRandomSeed () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This subroutine sets random seed for the pseudo-random number generator
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer :: i, n, clock
|
||||
integer, dimension(:), allocatable :: seed
|
||||
!-------------------------------------------------------------------------------------------
|
||||
call RANDOM_SEED ( size = n )
|
||||
allocate ( seed(n) )
|
||||
call SYSTEM_CLOCK ( COUNT = clock )
|
||||
seed = clock + 37 * (/ (i - 1, i = 1, n) /)
|
||||
call RANDOM_SEED ( PUT = seed )
|
||||
deallocate ( seed )
|
||||
end subroutine SetRandomSeed !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Generators for (random) properties of nanotubes
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8 function TubeMass ( i ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function returns the mass of the tube in Da
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4, intent(in) :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
TubeMass = M_2PI * RT(i) * LT(i) * TPBM * TPBD
|
||||
end function TubeMass !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function TubeSpecificHeat ( i ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function returns the specific heat of the tube
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4, intent(in) :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
TubeSpecificHeat = TPBSH
|
||||
end function TubeSpecificHeat !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Reading and printing of input parameters
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine LoadGoverningParameters () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function reads governing parameters from xdt file
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: Fuid, i
|
||||
character*512 :: Msg
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Fuid = OpenFile ( 'TMDGen.xdt', 'rt', '' )
|
||||
read ( unit = Fuid, fmt = '(e22.12)' ) LS0
|
||||
read ( unit = Fuid, fmt = '(e22.12)' ) HS0
|
||||
read ( unit = Fuid, fmt = '(e22.12)' ) DS0
|
||||
read ( unit = Fuid, fmt = '(i22)' ) BC_X0
|
||||
read ( unit = Fuid, fmt = '(i22)' ) BC_Y0
|
||||
read ( unit = Fuid, fmt = '(i22)' ) BC_Z0
|
||||
read ( unit = Fuid, fmt = '(i22)' ) ChiIndM
|
||||
read ( unit = Fuid, fmt = '(i22)' ) ChiIndN
|
||||
read ( unit = Fuid, fmt = '(e22.12)' ) LT0
|
||||
read ( unit = Fuid, fmt = '(i22)' ) SegType
|
||||
read ( unit = Fuid, fmt = '(i22)' ) NSeg0
|
||||
read ( unit = Fuid, fmt = '(e22.12)' ) LSeg0
|
||||
read ( unit = Fuid, fmt = '(e22.12)' ) DeltaT
|
||||
read ( unit = Fuid, fmt = '(i22)' ) NAmax
|
||||
read ( unit = Fuid, fmt = '(e22.12)' ) GeomPrec
|
||||
call CloseFile ( Fuid )
|
||||
end subroutine LoadGoverningParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine PrintGoverningParameters () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function prints governing parameters to xlg file
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: Fuid, i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Fuid = OpenFile ( 'TMDGen.xlg', 'wt', '' )
|
||||
write ( unit = Fuid, fmt = '(e22.12,a)' ) LS0, ' : LS0, Angstrom'
|
||||
write ( unit = Fuid, fmt = '(e22.12,a)' ) HS0, ' : HS0, Angstrom'
|
||||
write ( unit = Fuid, fmt = '(e22.12,a)' ) DS0, ' : DS0, g/cm**3'
|
||||
write ( unit = Fuid, fmt = '(e22.12,a)' ) DS0, ' : SC0, 1/A**2'
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) BC_X0, ' : BC_X0'
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) BC_Y0, ' : BC_Y0'
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) BC_Z0, ' : BC_Z0'
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) ChiIndM, ' : ChiIndM'
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) ChiIndN, ' : ChiIndN'
|
||||
write ( unit = Fuid, fmt = '(e22.12,a)' ) LT0, ' : LT0, Angstrom'
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) SegType, ' : SegType'
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) NSeg0, ' : NSeg0'
|
||||
write ( unit = Fuid, fmt = '(e22.12,a)' ) LSeg0, ' : LSeg0, Angstrom'
|
||||
write ( unit = Fuid, fmt = '(e22.12,a)' ) DeltaT, ' : DeltaT'
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) NAmax, ' : NAmax'
|
||||
write ( unit = Fuid, fmt = '(e22.12,a)' ) GeomPrec, ' : GeomPrec'
|
||||
call CloseFile ( Fuid )
|
||||
end subroutine PrintGoverningParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Printing of sample parameters
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine PrintSampleParameters ( ParType ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function prints the most imprtant parameters of the sample.
|
||||
! In the code, it used twice to print parameters of the desired and really generated samples.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
character*(*), intent(in) :: ParType
|
||||
real*8 :: MP, M, V
|
||||
!-------------------------------------------------------------------------------------------
|
||||
print '(a,a,a)', '*** ', trim(ParType), ' properties of the sample'
|
||||
print '(a34,a,f15.4,a)', 'L', ' : ', LS0, ' A'
|
||||
print '(a34,a,f15.4,a)', 'H', ' : ', HS0, ' A'
|
||||
print '(a34,a,f15.4,a)', 'Density', ' : ', DS0, ' g/cm**3'
|
||||
print '(a34,a,e15.8,a)', 'Volume', ' : ', VS0, ' A*3'
|
||||
print '(a34,a,e15.8,a)', 'Mass', ' : ', MS0, ' Da'
|
||||
print '(a34,a,i10)', 'BC_X', ' : ', BC_X0
|
||||
print '(a34,a,i10)', 'BC_Y', ' : ', BC_Y0
|
||||
print '(a34,a,i10)', 'BC_Z', ' : ', BC_Z0
|
||||
end subroutine PrintSampleParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Initializing of basic geometrical parameters of the generated sample
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine InitSample () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function initializes the geometrical parameters of the sample (sizes, etc.)
|
||||
!-------------------------------------------------------------------------------------------
|
||||
|
||||
BC_X = BC_X0
|
||||
BC_Y = BC_Y0
|
||||
BC_Z = BC_Z0
|
||||
DomXmin = - LS0 / 2.0d+00
|
||||
DomXmax = LS0 / 2.0d+00
|
||||
DomYmin = - LS0 / 2.0d+00
|
||||
DomYmax = LS0 / 2.0d+00
|
||||
DomZmin = - HS0 / 2.0d+00
|
||||
DomZmax = HS0 / 2.0d+00
|
||||
|
||||
if ( BC_X0 == 0 ) then
|
||||
DomXmin = 0.0d+00
|
||||
DomXmax = LS0
|
||||
end if
|
||||
if ( BC_Y0 == 0 ) then
|
||||
DomYmin = 0.0d+00
|
||||
DomYmax = LS0
|
||||
end if
|
||||
if ( BC_Z0 == 0 ) then
|
||||
DomZmin = 0.0d+00
|
||||
DomZmax = HS0
|
||||
end if
|
||||
|
||||
DomLX = DomXmax - DomXmin
|
||||
DomLY = DomYmax - DomYmin
|
||||
DomLZ = DomZmax - DomZmin
|
||||
DomLXHalf = 0.5d+00 * DomLX
|
||||
DomLYHalf = 0.5d+00 * DomLY
|
||||
DomLZHalf = 0.5d+00 * DomLZ
|
||||
|
||||
DS0 = DS0 / ( K_MDDU / 1.0d+03 )
|
||||
VS0 = LS0 * LS0 * HS0
|
||||
MS0 = DS0 * VS0
|
||||
end subroutine InitSample !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! A few auxiliary functions
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine GetTubeEnds ( X0, X1, i ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function calculates coordinates of two ends of nanotube i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2), intent(out) :: X0, X1
|
||||
integer*4, intent(in) :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: LT2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
LT2 = 0.5d+00 * LT(i)
|
||||
X0 = CT(i,0:2) - LT2 * DT(i,0:2)
|
||||
X1 = CT(i,0:2) + LT2 * DT(i,0:2)
|
||||
end subroutine GetTubeEnds !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
logical function IsTubeInside ( i ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function returns true if nanotube i lies inside the sample. Otherwise it returns false.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4, intent(in) :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: n
|
||||
real*8, dimension(0:2) :: X0, X1, Xmin, Xmax
|
||||
!-------------------------------------------------------------------------------------------
|
||||
IsTubeInside = .true.
|
||||
if ( BC_X == 1 .and. BC_Y == 1 .and. BC_Z == 1 ) return
|
||||
call GetTubeEnds ( X0, X1, i )
|
||||
do n = 0, 2
|
||||
Xmin(n) = min ( X0(n), X1(n) )
|
||||
Xmax(n) = max ( X0(n), X1(n) )
|
||||
end do
|
||||
IsTubeInside = .false.
|
||||
if ( BC_X == 0 .and. ( Xmin(0) < DomXmin .or. Xmax(0) > DomXmax ) ) return
|
||||
if ( BC_Y == 0 .and. ( Xmin(1) < DomYmin .or. Xmax(1) > DomYmax ) ) return
|
||||
if ( BC_Z == 0 .and. ( Xmin(2) < DomZmin .or. Xmax(2) > DomZmax ) ) return
|
||||
IsTubeInside = .true.
|
||||
end function IsTubeInside !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module TMDGenData !*****************************************************************************
|
||||
@ -1,45 +0,0 @@
|
||||
newton on
|
||||
log cnt.log
|
||||
echo both
|
||||
units metal
|
||||
lattice sc 1.0
|
||||
boundary p p fs
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 5 delay 0 check yes
|
||||
|
||||
atom_style cnt
|
||||
#cut, RT, STRMode, BendingMode, STRParams, YMType, TPMType, TPMSSTP.xrs, TPMA.xrs
|
||||
pair_style cnt/cnt 45.0 6.785 1 0 3 0 0 ../../../potentials/TPMSSTP.xrs ../../../potentials/TPMA.xrs
|
||||
read_data TMDSample.init
|
||||
pair_coeff * *
|
||||
|
||||
velocity all create 600.0 2019
|
||||
timestep 0.010
|
||||
fix 1 all nve
|
||||
#fix 1 all nvt temp 300.0 300.0 1.0
|
||||
thermo_modify flush yes
|
||||
thermo 1
|
||||
reset_timestep 0
|
||||
|
||||
compute Es all cnt/Es
|
||||
compute Eb all cnt/Eb
|
||||
compute Et all cnt/Et
|
||||
compute Ek all ke/atom
|
||||
compute Es_tot all cnt/Es_tot
|
||||
compute Eb_tot all cnt/Eb_tot
|
||||
compute Et_tot all cnt/Et_tot
|
||||
compute Ep_tot all pe
|
||||
compute Ek_tot all ke
|
||||
variable time_ equal time
|
||||
variable Ep_ equal c_Ep_tot
|
||||
variable Ek_ equal c_Ek_tot
|
||||
variable Etot_ equal v_Ek_+v_Ep_
|
||||
variable Es_ equal c_Es_tot
|
||||
variable Eb_ equal c_Eb_tot
|
||||
variable Et_ equal c_Et_tot
|
||||
|
||||
dump out_dump all custom 50 config_E.dump id type x y z c_Es c_Eb c_Et c_Ek ix iy iz
|
||||
fix out_info all print 10 "${time_} ${Etot_} ${Ek_} ${Ep_} ${Es_} ${Eb_} ${Et_}" file "E.txt" screen no
|
||||
|
||||
run 50
|
||||
write_data system_E.data
|
||||
@ -1,144 +0,0 @@
|
||||
module TPMGeom !************************************************************************************
|
||||
!
|
||||
! Geometry functions for TPM force field
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, 2020, Version 13.00
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
use TPMLib
|
||||
|
||||
implicit none
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Constants
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4, parameter :: MD_LINES_NONPAR = 0
|
||||
integer*4, parameter :: MD_LINES_PAR = 1
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Global variables
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
! Coordinates of the whole domain
|
||||
real*8 :: DomXmin, DomXmax, DomYmin, DomYmax, DomZmin, DomZmax
|
||||
real*8 :: DomLX, DomLY, DomLZ
|
||||
real*8 :: DomLXhalf, DomLYhalf, DomLZhalf
|
||||
|
||||
! Boundary conditions
|
||||
integer*4 :: BC_X = 0
|
||||
integer*4 :: BC_Y = 0
|
||||
integer*4 :: BC_Z = 0
|
||||
|
||||
! Skin parameter in NBL and related algorithms
|
||||
real*8 :: Rskin = 1.0d+00
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
subroutine ApplyPeriodicBC ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This subroutine changes coortinates of the point accorning to periodic boundary conditions
|
||||
! it order to makesure that the point is inside the computational cell
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2), intent(inout) :: R
|
||||
!-------------------------------------------------------------------------------------------
|
||||
! These commented lines implemment the more general, but less efficient algorithm
|
||||
!if ( BC_X == 1 ) R(0) = R(0) - DomLX * roundint ( R(0) / DomLX )
|
||||
!if ( BC_Y == 1 ) R(1) = R(1) - DomLY * roundint ( R(1) / DomLY )
|
||||
!if ( BC_Z == 1 ) R(2) = R(2) - DomLZ * roundint ( R(2) / DomLZ )
|
||||
if ( BC_X == 1 ) then
|
||||
if ( R(0) .GT. DomLXHalf ) then
|
||||
R(0) = R(0) - DomLX
|
||||
else if ( R(0) .LT. - DomLXHalf ) then
|
||||
R(0) = R(0) + DomLX
|
||||
end if
|
||||
end if
|
||||
if ( BC_Y == 1 ) then
|
||||
if ( R(1) .GT. DomLYHalf ) then
|
||||
R(1) = R(1) - DomLY
|
||||
else if ( R(1) .LT. - DomLYHalf ) then
|
||||
R(1) = R(1) + DomLY
|
||||
end if
|
||||
end if
|
||||
if ( BC_Z == 1 ) then
|
||||
if ( R(2) .GT. DomLZHalf ) then
|
||||
R(2) = R(2) - DomLZ
|
||||
else if ( R(2) .LT. - DomLZHalf ) then
|
||||
R(2) = R(2) + DomLZ
|
||||
end if
|
||||
end if
|
||||
end subroutine ApplyPeriodicBC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine LinePoint ( Displacement, Q, R1, L1, R0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function calculates the point Q of projection of point R0 on line (R1,L1)
|
||||
! Q = R1 + Disaplacement * L1
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(inout) :: Displacement
|
||||
real*8, dimension(0:2), intent(inout) :: Q
|
||||
real*8, dimension(0:2), intent(in) :: R1, L1, R0
|
||||
!--------------------------------------------------------------------------------------------
|
||||
Q = R0 - R1
|
||||
! Here we take into account periodic boundaries
|
||||
call ApplyPeriodicBC ( Q )
|
||||
Displacement = S_V3xV3 ( Q, L1 )
|
||||
Q = R1 + Displacement * L1
|
||||
end subroutine LinePoint !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
integer*4 function LineLine ( H, cosA, D1, D2, L12, R1, L1, R2, L2, Prec ) !!!!!!!!!!!!!!!!!
|
||||
! This function determines the neares distance H between two lines (R1,L1) and (R2,L2)
|
||||
!-------------------------------------------------------------------------------------------
|
||||
! Input values:
|
||||
! R1, L1, point and direction of line 1
|
||||
! R2, L2, point and direction of line 2
|
||||
! Prec, precision for the case L1 * L2 = 0 (parallel lines)
|
||||
! Return values:
|
||||
! H, minimal distance between lines
|
||||
! cosA, cosine of angle between lines
|
||||
! D1, D2, displacemets
|
||||
! L12, unit vector directed along the closes distance
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(inout) :: H, cosA, D1, D2
|
||||
real*8, dimension(0:2), intent(out) :: L12
|
||||
real*8, dimension(0:2), intent(in) :: R1, L1, R2, L2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(in) :: Prec
|
||||
real*8, dimension(0:2) :: Q1, Q2, R
|
||||
real*8 :: C, DD1, DD2, C1, C2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
cosA = S_V3xV3 ( L1, L2 )
|
||||
C = 1.0 - sqr ( cosA )
|
||||
if ( C < Prec ) then ! Lines are parallel to each other
|
||||
LineLine = MD_LINES_PAR
|
||||
return
|
||||
end if
|
||||
LineLine = MD_LINES_NONPAR
|
||||
R = R2 - R1
|
||||
! Here we take into account periodic boundaries
|
||||
call ApplyPeriodicBC ( R )
|
||||
DD1 = S_V3xV3 ( R, L1 )
|
||||
DD2 = S_V3xV3 ( R, L2 )
|
||||
D1 = ( cosA * DD2 - DD1 ) / C
|
||||
D2 = ( DD2 - cosA * DD1 ) / C
|
||||
Q1 = R1 - D1 * L1
|
||||
Q2 = R2 - D2 * L2
|
||||
L12 = Q2 - Q1
|
||||
! Here we take into account periodic boundaries
|
||||
call ApplyPeriodicBC ( L12 )
|
||||
H = S_V3norm3 ( L12 )
|
||||
if ( H < Prec ) then ! Lines intersect each other
|
||||
C1 = signum ( D1 )
|
||||
C2 = signum ( D1 ) * signum ( cosA )
|
||||
Q1 = C1 * L1
|
||||
Q2 = C2 * L2
|
||||
call V3_V3xxV3 ( L12, Q1, Q2 )
|
||||
call V3_ort ( L12 )
|
||||
else ! No intersection
|
||||
L12 = L12 / H
|
||||
end if
|
||||
end function LineLine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module TPMGeom !********************************************************************************
|
||||
@ -1,205 +0,0 @@
|
||||
module TPMLib !*************************************************************************************
|
||||
!
|
||||
! Common constants, types, and functions for TPM force field
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, 2020, Version 13.00
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
implicit none
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Mathematical constants
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8, parameter :: M_PI_2 = 1.57079632679489661923
|
||||
real*8, parameter :: M_PI = 3.14159265358979323846
|
||||
real*8, parameter :: M_3PI_2 = 4.71238898038468985769
|
||||
real*8, parameter :: M_2PI = 6.28318530717958647692
|
||||
real*8, parameter :: M_PI_180 = 0.017453292519943295769
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Physical unit constants
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8, parameter :: K_AMU = 1.66056E-27 ! a.m.u. (atomic mass unit, Dalton)
|
||||
real*8, parameter :: K_EV = 1.60217646e-19 ! eV (electron-volt)
|
||||
|
||||
real*8, parameter :: K_MDLU = 1.0E-10 ! MD length unit (m)
|
||||
real*8, parameter :: K_MDEU = K_EV ! MD energy unit (J)
|
||||
real*8, parameter :: K_MDMU = K_AMU ! MD mass unit (kg)
|
||||
real*8, parameter :: K_MDFU = K_MDEU / K_MDLU ! MD force unit (N)
|
||||
real*8, parameter :: K_MDCU = K_MDEU / K_MDMU ! MD specific heat unit (J/(kg*K))
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Global variables
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4 :: StdUID = 31
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Simple mathematical functions
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8 function rad ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: X
|
||||
!-------------------------------------------------------------------------------------------
|
||||
rad = X * M_PI_180
|
||||
end function rad !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function sqr ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: X
|
||||
!-------------------------------------------------------------------------------------------
|
||||
sqr = X * X
|
||||
end function sqr !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
integer*4 function signum ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: X
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( X > 0 ) then
|
||||
signum = 1
|
||||
else if ( X < 0 ) then
|
||||
signum = -1
|
||||
else
|
||||
signum = 0
|
||||
end if
|
||||
end function signum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Vector & matrix functions
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8 function S_V3xx ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, dimension(0:2), intent(in) :: V
|
||||
!-------------------------------------------------------------------------------------------
|
||||
S_V3xx = V(0) * V(0) + V(1) * V(1) + V(2) * V(2)
|
||||
end function S_V3xx !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function S_V3xV3 ( V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, dimension(0:2), intent(in) :: V1, V2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
S_V3xV3 = V1(0) * V2(0) + V1(1) * V2(1) + V1(2) * V2(2)
|
||||
end function S_V3xV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function S_V3norm3 ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, dimension(0:2), intent(in) :: V
|
||||
!-------------------------------------------------------------------------------------------
|
||||
S_V3norm3 = dsqrt ( V(0) * V(0) + V(1) * V(1) + V(2) * V(2) )
|
||||
end function S_V3norm3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine V3_ort ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! Vector production
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2), intent(inout) :: V
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: Vabs
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Vabs = S_V3norm3 ( V )
|
||||
V(0) = V(0) / Vabs
|
||||
V(1) = V(1) / Vabs
|
||||
V(2) = V(2) / Vabs
|
||||
end subroutine V3_ort !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine V3_V3xxV3 ( V, V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! Vector production
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2), intent(out) :: V
|
||||
real*8, dimension(0:2), intent(in) :: V1, V2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
V(0) = V1(1) * V2(2) - V1(2) * V2(1)
|
||||
V(1) = V1(2) * V2(0) - V1(0) * V2(2)
|
||||
V(2) = V1(0) * V2(1) - V1(1) * V2(0)
|
||||
end subroutine V3_V3xxV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Handling of spherical and Euler angles
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine RotationMatrix3 ( M, Psi, Tet, Phi ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! Ksi, Tet and Phi are Euler angles
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2,0:2), intent(out) :: M
|
||||
real*8, intent(in) :: Psi, Tet, Phi
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: cK, sK, cT, sT, cP, sP
|
||||
!-------------------------------------------------------------------------------------------
|
||||
cK = dcos ( Psi )
|
||||
sK = dsin ( Psi )
|
||||
cT = dcos ( Tet )
|
||||
sT = dsin ( Tet )
|
||||
cP = dcos ( Phi )
|
||||
sP = dsin ( Phi )
|
||||
M(0,0) = cP * cK - sK * sP * cT
|
||||
M(0,1) = cP * sK + sP * cT * cK
|
||||
M(0,2) = sP * sT
|
||||
M(1,0) = - sP * cK - cP * cT * sK
|
||||
M(1,1) = - sP * sK + cP * cT * cK
|
||||
M(1,2) = cP * sT
|
||||
M(2,0) = sT * sK
|
||||
M(2,1) = - sT * cK
|
||||
M(2,2) = cT
|
||||
end subroutine RotationMatrix3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine EulerAngles ( Psi, Tet, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(out) :: Tet, Psi
|
||||
real*8, dimension(0:2), intent(in) :: L
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Tet = acos ( L(2) )
|
||||
Psi = atan2 ( L(1), L(0) )
|
||||
if ( Psi > M_3PI_2 ) then
|
||||
Psi = Psi - M_3PI_2
|
||||
else
|
||||
Psi = Psi + M_PI_2
|
||||
end if
|
||||
end subroutine EulerAngles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! File inout and output
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4 function OpenFile ( Name, Params, Path ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
character*(*), intent(in) :: Name, Params, Path
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: Fuid
|
||||
character*512 :: FullName, Msg, Name1, Action1, Status1, Form1, Position1
|
||||
!-------------------------------------------------------------------------------------------
|
||||
OpenFile = StdUID + 1
|
||||
if ( Params(1:1) == 'r' ) then
|
||||
Action1 = 'read'
|
||||
Status1 = 'old'
|
||||
Position1 = 'rewind'
|
||||
else if ( Params(1:1) == 'w' ) then
|
||||
Action1 = 'write'
|
||||
Status1 = 'replace'
|
||||
Position1 = 'rewind'
|
||||
else if ( Params(1:1) == 'a' ) then
|
||||
Action1 = 'write'
|
||||
Status1 = 'old'
|
||||
Position1 = 'append'
|
||||
endif
|
||||
if ( Params(2:2) == 'b' ) then
|
||||
Form1 = 'binary'
|
||||
else
|
||||
Form1 = 'formatted'
|
||||
endif
|
||||
open ( unit = OpenFile, file = Name, form = Form1, action = Action1, status = Status1, position = Position1 )
|
||||
StdUID = StdUID + 1
|
||||
return
|
||||
end function OpenFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine CloseFile ( Fuid ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
integer*4, intent(inout) :: Fuid
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( Fuid < 0 ) return
|
||||
close ( unit = Fuid )
|
||||
Fuid = -1
|
||||
end subroutine CloseFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module TPMLib !*********************************************************************************
|
||||
|
||||
@ -1,98 +0,0 @@
|
||||
module LinFun2 !************************************************************************************
|
||||
!
|
||||
! Bi-linear functions and their derivatives.
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
implicit none
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
real*8 function CalcLinFun1_0 ( i, X, N, P, F ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
integer*4, intent(in) :: i, N
|
||||
real*8, intent(in) :: X
|
||||
real*8, dimension(0:N-1), intent(in) :: P
|
||||
real*8, dimension(0:N-1), intent(inout) :: F
|
||||
integer*4 :: i1
|
||||
real*8 :: A, A0
|
||||
!-------------------------------------------------------------------------------------------
|
||||
i1 = i - 1
|
||||
A0 = ( P(i) - X ) / ( P(i) - P(i1) )
|
||||
A = 1.0d+00 - A0
|
||||
CalcLinFun1_0 = A0 * F(i1) + A * F(i)
|
||||
end function CalcLinFun1_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine CalcLinFun1_1 ( S, Sx1, i, X, N, P, F, Fx ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(out) :: S, Sx1
|
||||
integer*4, intent(in) :: i, N
|
||||
real*8, intent(in) :: X
|
||||
real*8, dimension(0:N-1), intent(in) :: P
|
||||
real*8, dimension(0:N-1), intent(inout) :: F, Fx
|
||||
integer*4 :: i1
|
||||
real*8 :: A, A0
|
||||
!-------------------------------------------------------------------------------------------
|
||||
i1 = i - 1
|
||||
A0 = ( P(i) - X ) / ( P(i) - P(i1) )
|
||||
A = 1.0d+00 - A0
|
||||
S = A0 * F(i1) + A * F(i)
|
||||
Sx1 = A0 * Fx(i1) + A * Fx(i)
|
||||
end subroutine CalcLinFun1_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function CalcLinFun2_0 ( i, j, X, Y, N1, N2, P1, P2, F ) !!
|
||||
integer*4, intent(in) :: i, j, N1, N2
|
||||
real*8, intent(in) :: X, Y
|
||||
real*8, dimension(0:N1-1), intent(in) :: P1
|
||||
real*8, dimension(0:N2-1), intent(in) :: P2
|
||||
real*8, dimension(0:N1-1,0:N2-1), intent(inout) :: F
|
||||
integer*4 :: i1, j1
|
||||
real*8 :: A, A0, B, B0, G, G0
|
||||
!-------------------------------------------------------------------------------------------
|
||||
i1 = i - 1
|
||||
j1 = j - 1
|
||||
A0 = ( P1(i) - X ) / ( P1(i) - P1(i1) )
|
||||
A = 1.0d+00 - A0
|
||||
B0 = ( P2(j) - Y ) / ( P2(j) - P2(j1) )
|
||||
B = 1.0d+00 - B0
|
||||
G = B0 * F(i,j1) + B * F(i,j)
|
||||
G0 = B0 * F(i1,j1) + B * F(i1,j)
|
||||
CalcLinFun2_0 = A0 * G0 + A * G
|
||||
end function CalcLinFun2_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine CalcLinFun2_1 ( S, Sx1, Sy1, i, j, X, Y, N1, N2, P1, P2, F, Fx, Fy ) !!!!!!!!!!!!
|
||||
real*8, intent(out) :: S, Sx1, Sy1
|
||||
integer*4, intent(in) :: i, j, N1, N2
|
||||
real*8, intent(in) :: X, Y
|
||||
real*8, dimension(0:N1-1), intent(in) :: P1
|
||||
real*8, dimension(0:N2-1), intent(in) :: P2
|
||||
real*8, dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fx, Fy
|
||||
integer*4 :: i1, j1
|
||||
real*8 :: A, A0, B, B0, G, G0
|
||||
!-------------------------------------------------------------------------------------------
|
||||
i1 = i - 1
|
||||
j1 = j - 1
|
||||
A0 = ( P1(i) - X ) / ( P1(i) - P1(i1) )
|
||||
A = 1.0d+00 - A0
|
||||
B0 = ( P2(j) - Y ) / ( P2(j) - P2(j1) )
|
||||
B = 1.0d+00 - B0
|
||||
|
||||
G = B0 * F(i,j1) + B * F(i,j)
|
||||
G0 = B0 * F(i1,j1) + B * F(i1,j)
|
||||
S = A0 * G0 + A * G
|
||||
|
||||
G = B0 * Fx(i,j1) + B * Fx(i,j)
|
||||
G0 = B0 * Fx(i1,j1) + B * Fx(i1,j)
|
||||
Sx1 = A0 * G0 + A * G
|
||||
|
||||
G = B0 * Fy(i,j1) + B * Fy(i,j)
|
||||
G0 = B0 * Fy(i1,j1) + B * Fy(i1,j)
|
||||
Sy1 = A0 * G0 + A * G
|
||||
|
||||
end subroutine CalcLinFun2_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module LinFun2 !********************************************************************************
|
||||
@ -1,35 +0,0 @@
|
||||
#---------------------------------------------------------------------------------------------------
|
||||
#
|
||||
# This is Makefile for builing the executable TMDPotGen
|
||||
#
|
||||
# Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, 2020, Version 13.00
|
||||
#
|
||||
#---------------------------------------------------------------------------------------------------
|
||||
EXEPATH = .
|
||||
|
||||
F90 = ifort
|
||||
F90FLAGS = -Ofast -mcmodel=medium
|
||||
#F90 = pgf90
|
||||
#F90FLAGS = -fast -mcmodel=medium
|
||||
LDFLAGS =
|
||||
|
||||
OBJS = TPMLib.o LinFun2.o Spline1.o Spline2.o TPMGeom.o TubePotBase.o TubePotTrue.o TubePotMono.o TMDPotGen.o
|
||||
EXE = $(EXEPATH)/TMDPotGen
|
||||
|
||||
# compile and load
|
||||
default:
|
||||
@echo " "
|
||||
@echo "Compiling Code of Program TMDPotGen"
|
||||
@echo "FORTRAN 90"
|
||||
$(MAKE) $(EXE)
|
||||
|
||||
$(EXE): $(OBJS)
|
||||
$(F90) $(F90FLAGS) $(LDFLAGS) -o $(EXE) $(OBJS)
|
||||
|
||||
.SUFFIXES: .f90 .o
|
||||
.f90.o:
|
||||
$(F90) $(F90FLAGS) -c $*.f90
|
||||
|
||||
clean:
|
||||
rm -f *.o
|
||||
|
||||
@ -1,177 +0,0 @@
|
||||
module Spline1 !************************************************************************************
|
||||
!
|
||||
! One-dimensional cubic spline function.
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
implicit none
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
real*8 function ValueSpline1_0 ( X, Xi, Xi_1, Yi, Yi_1, Mi, Mi_1, Hi_1 ) !!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: X, Xi, Xi_1, Yi, Yi_1, Mi, Mi_1, Hi_1
|
||||
real*8 :: H26, HL, HR
|
||||
!-------------------------------------------------------------------------------------------
|
||||
H26 = Hi_1 * Hi_1 / 6.0
|
||||
Hl = X - Xi_1
|
||||
Hr = Xi - X
|
||||
ValueSpline1_0 = ( ( Mi_1 * Hr * Hr * Hr + Mi * Hl * Hl * Hl ) / 6.0 + ( Yi_1 - Mi_1 * H26 ) * Hr &
|
||||
+ ( Yi - Mi * H26 ) * Hl ) / Hi_1
|
||||
end function ValueSpline1_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine ValueSpline1_1 ( S, S1, X, Xi, Xi_1, Yi, Yi_1, Mi, Mi_1, Hi_1 ) !!!!!!!!!!!!!!!!!
|
||||
real*8, intent(out) :: S, S1
|
||||
real*8, intent(in) :: X, Xi, Xi_1, Yi, Yi_1, Mi, Mi_1, Hi_1
|
||||
real*8 :: H6, H26, HL, HR, HL2, HR2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
H6 = Hi_1 / 6.0d+00
|
||||
H26 = Hi_1 * H6
|
||||
HL = X - Xi_1
|
||||
HR = Xi - X
|
||||
HL2 = HL * HL
|
||||
HR2 = HR * HR
|
||||
S = ( ( Mi_1 * HR2 * Hr + Mi * HL2 * Hl ) / 6.0 + ( Yi_1 - Mi_1 * H26 ) * HR + ( Yi - Mi * H26 ) * HL ) / Hi_1
|
||||
S1 = ( ( Mi * HL2 - Mi_1 * HR2 ) / 2.0d+00 + Yi - Yi_1 ) / Hi_1 - H6 * ( Mi - Mi_1 )
|
||||
end subroutine ValueSpline1_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine sprogonka3 ( N, K0, K1, K2, F, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
integer*4, intent(in) :: N
|
||||
real*8, dimension(0:N-1), intent(in) :: K0, K1, K2
|
||||
real*8, dimension(0:N-1), intent(inout) :: F, X
|
||||
real*8 :: D
|
||||
integer*4 :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
X(0) = F(0) / K1(0)
|
||||
F(0) = - K2(0) / K1(0)
|
||||
do i = 1, N - 1
|
||||
D = - ( K1(i) + F(i-1) * K0(i) )
|
||||
X(i) = ( K0(i) * X(i-1) - F(i) ) / D
|
||||
F(i) = K2(i) / D
|
||||
end do
|
||||
do i = N - 2, 0, -1
|
||||
X(i) = X(i) + F(i) * X(i+1)
|
||||
end do
|
||||
end subroutine sprogonka3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine CreateSpline1 ( CL, CR, N, P, F, M, D, K0, K1, K2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
integer*4, intent(in) :: CL, CR, N
|
||||
real*8, dimension (0:N-1), intent(in) :: P, F
|
||||
real*8, dimension (0:N-1), intent(inout):: M, D, K0, K1, K2
|
||||
integer*4 :: i
|
||||
real*8 :: Z
|
||||
!-------------------------------------------------------------------------------------------
|
||||
do i = 1, N - 1
|
||||
K0(i) = P(i) - P(i-1)
|
||||
K1(i) = ( F(i) - F(i-1) ) / K0(i)
|
||||
end do
|
||||
select case ( CL )
|
||||
case (1)
|
||||
K1(0) = 2.0d+00 / 3.0d+00
|
||||
K2(0) = 1.0d+00 / 3.0d+00
|
||||
D (0) = 2 * ( K1(1) - M(0) ) / K0(1)
|
||||
case (2)
|
||||
K1(0) = 1.0d+00
|
||||
K2(0) = 0.0d+00
|
||||
D(0) = M(0)
|
||||
case (3)
|
||||
K1(0) = 1.0d+00
|
||||
K2(0) = 0.0d+00
|
||||
D(0) = 0.0d+00
|
||||
end select
|
||||
Z = K1(N-1)
|
||||
do i = 1, N - 2
|
||||
D(i) = 6.0d+00 * ( K1(i+1) - K1(i) )
|
||||
K2(i) = K0(i+1)
|
||||
K1(i) = 2.0d+00 * ( K2(i) + K0(i) )
|
||||
end do
|
||||
select case ( CR )
|
||||
case (1)
|
||||
D(N-1) = 2.0d+00 * ( M(N-1) - Z ) / K0(N-1)
|
||||
K1(N-1) = 2.0d+00 / 3.0d+00
|
||||
K0(N-1) = 1.0d+00 / 3.0d+00
|
||||
case (2)
|
||||
K1(N-1) = 1.0d+00
|
||||
K0(N-1) = 0.0d+00
|
||||
D(N-1) = M(N-1)
|
||||
case (3)
|
||||
K1(N-1) = 1.0d+00
|
||||
K0(N-1) = 0.0d+00
|
||||
D(N-1) = 0.0d+00
|
||||
end select
|
||||
call sprogonka3 ( N, K0, K1, K2, D, M )
|
||||
end subroutine CreateSpline1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function CalcSpline1_0 ( i, X, N, P, F, M ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
integer*4, intent(in) :: i, N
|
||||
real*8, intent(in) :: X
|
||||
real*8, dimension(0:N-1), intent(in) :: P, F, M
|
||||
integer*4 :: j
|
||||
real*8 :: HL, HR, H, H6, H26, HR2, HL2, HRH, HLH
|
||||
!-------------------------------------------------------------------------------------------
|
||||
j = i - 1
|
||||
HL = X - P(j)
|
||||
HR = P(i) - X
|
||||
H = P(i) - P(j)
|
||||
H6 = H / 6.0d+00
|
||||
H26 = H * H6
|
||||
HL2 = HL * HL
|
||||
HR2 = HR * HR
|
||||
HLH = HL / H
|
||||
HRH = HR / H
|
||||
CalcSpline1_0 = ( M(j) * HR2 * HRH + M(i) * HL2 * HLH ) / 6.0d+00 + ( F(j) - M(j) * H26 ) * HRH &
|
||||
+ ( F(i) - M(i) * H26 ) * HLH
|
||||
end function CalcSpline1_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine CalcSpline1_1 ( S, S1, i, X, N, P, F, M ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(out) :: S, S1
|
||||
integer*4, intent(in) :: i, N
|
||||
real*8, intent(in) :: X
|
||||
real*8, dimension(0:N-1), intent(in) :: P, F, M
|
||||
integer*4 :: j
|
||||
real*8 :: HL, HR, H, H6, H26, HR2, HL2, HRH, HLH
|
||||
!-------------------------------------------------------------------------------------------
|
||||
j = i - 1
|
||||
HL = X - P(j)
|
||||
HR = P(i) - X
|
||||
H = P(i) - P(j)
|
||||
H6 = H / 6.0d+00
|
||||
H26 = H * H6
|
||||
HL2 = HL * HL
|
||||
HR2 = HR * HR
|
||||
HLH = HL / H
|
||||
HRH = HR / H
|
||||
S = ( M(j) * HR2 * HRH + M(i) * HL2 * HLH ) / 6.0d+00 + ( F(j) - M(j) * H26 ) * HRH + ( F(i) - M(i) * H26 ) * HLH
|
||||
S1 = ( ( M(i) * HL2 - M(j) * HR2 ) / 2.0d+00 + F(i) - F(j) ) / H - H6 * ( M(i) - M(j) )
|
||||
end subroutine CalcSpline1_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine CalcSpline1_2 ( S, S1, S2, i, X, N, P, F, M ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(out) :: S, S1, S2
|
||||
integer*4, intent(in) :: i, N
|
||||
real*8, intent(in) :: X
|
||||
real*8, dimension(0:N-1), intent(in) :: P, F, M
|
||||
integer*4 :: j
|
||||
real*8 :: HL, HR, H, H6, H26, HR2, HL2, HRH, HLH
|
||||
!-------------------------------------------------------------------------------------------
|
||||
j = i - 1
|
||||
HL = X - P(j)
|
||||
HR = P(i) - X
|
||||
H = P(i) - P(j)
|
||||
H6 = H / 6.0d+00
|
||||
H26 = H * H6
|
||||
HL2 = HL * HL
|
||||
HR2 = HR * HR
|
||||
HLH = HL / H
|
||||
HRH = HR / H
|
||||
S = ( M(j) * HR2 * HRH + M(i) * HL2 * HLH ) / 6.0d+00 + ( F(j) - M(j) * H26 ) * HRH + ( F(i) - M(i) * H26 ) * HLH
|
||||
S1 = ( ( M(i) * HL2 - M(j) * HR2 ) / 2.0d+00 + F(i) - F(j) ) / H - H6 * ( M(i) - M(j) )
|
||||
S2 = M(j) * HRH + M(i) * HLH
|
||||
end subroutine CalcSpline1_2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module Spline1 !********************************************************************************
|
||||
|
||||
@ -1,171 +0,0 @@
|
||||
module Spline2 !************************************************************************************
|
||||
!
|
||||
! Two-dimensional cubic spline function.
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
use Spline1
|
||||
|
||||
implicit none
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
subroutine CreateSpline2 ( CL, CD, CR, CU, N1, N2, N, P1, P2, F, Fxx, Fyy, Fxxyy, FF, MM, DD, K0, K1, K2 )
|
||||
integer*4, intent(in) :: CL, CD, CR, CU, N1, N2, N
|
||||
real*8, dimension(0:N1-1), intent(in) :: P1
|
||||
real*8, dimension(0:N2-1), intent(in) :: P2
|
||||
real*8, dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy
|
||||
real*8, dimension(0:N-1), intent(inout) :: FF, MM, DD, K0, K1, K2
|
||||
integer*4 :: II
|
||||
!-------------------------------------------------------------------------------------------
|
||||
do II = 0, N2 - 1
|
||||
FF(0:N1-1) = F(0:N1-1,II)
|
||||
MM(0) = Fxx(0,II)
|
||||
MM(N1-1) = Fxx(N1-1,II)
|
||||
call CreateSpline1 ( CL, CR, N1, P1, FF, MM, DD, K0, K1, K2 )
|
||||
Fxx(0:N1-1,II) = MM(0:N1-1)
|
||||
end do
|
||||
do II = 0, N1 - 1
|
||||
MM(0) = Fyy(II,0)
|
||||
MM(N-1) = Fyy(II,N2-1)
|
||||
FF(0:N2-1) = F(II,0:N2-1)
|
||||
call CreateSpline1 ( CD, CU, N2, P2, FF, MM, DD, K0, K1, K2 )
|
||||
Fyy(II,0:N2-1) = MM(0:N2-1)
|
||||
end do
|
||||
FF(0:N1-1) = Fyy(0:N1-1,0 )
|
||||
call CreateSpline1 ( 3, 3, N1, P1, FF, MM, DD, K0, K1, K2 )
|
||||
Fxxyy(0:N1-1,0) = MM(0:N1-1)
|
||||
FF(0:N1-1) = Fyy(0:N1-1,N2-1 )
|
||||
call CreateSpline1 ( 3, 3, N1, P1, FF, MM, DD, K0, K1, k2 )
|
||||
Fxxyy(0:N1-1,N2-1) = MM(0:N1-1)
|
||||
do II = 1, N1 - 2
|
||||
MM(0) = Fxxyy(II,0)
|
||||
MM(N-1) = Fxxyy(II,N2-1)
|
||||
FF(0:N2-1) = Fxx(II,0:N2-1)
|
||||
call CreateSpline1 ( 2 , 2, N2, P2, FF, MM, DD, K0, K1, K2 )
|
||||
Fxxyy(II,0:N2-1) = MM(0:N2-1)
|
||||
end do
|
||||
end subroutine CreateSpline2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine CreateSpline2Ext ( CL, CD, CR, CU, N1, N1A, N2, N2A, N, P1, P2, F, Fxx, Fyy, Fxxyy, FF, MM, DD, K0, K1, K2 )
|
||||
integer*4, intent(in) :: CL, CD, CR, CU, N1, N1A, N2, N2A, N
|
||||
real*8, dimension(0:N1-1), intent(in) :: P1
|
||||
real*8, dimension(0:N2-1), intent(in) :: P2
|
||||
real*8, dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy
|
||||
real*8, dimension(0:N-1), intent(inout) :: FF, MM, DD, K0, K1, K2
|
||||
integer*4 :: II
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Fxx = 0.0d+00
|
||||
Fyy = 0.0d+00
|
||||
Fxxyy = 0.0d+00
|
||||
|
||||
do II = 0, N2A
|
||||
FF(0:N1-1) = F(0:N1-1,II)
|
||||
MM(0) = Fxx(0,II)
|
||||
MM(N1-1) = Fxx(N1-1,II)
|
||||
call CreateSpline1 ( CL, CR, N1, P1, FF, MM, DD, K0, K1, K2 )
|
||||
Fxx(0:N1-1,II) = MM(0:N1-1)
|
||||
end do
|
||||
|
||||
do II = N2A + 1, N2 - 1
|
||||
FF(0:N1-N1A-1) = F(N1A:N1-1,II)
|
||||
MM(0) = Fxx(N1A,II)
|
||||
MM(N1-N1A-1) = Fxx(N1-1,II)
|
||||
call CreateSpline1 ( CL, CR, N1 - N1A, P1, FF, MM, DD, K0, K1, K2 )
|
||||
Fxx(N1A:N1-1,II) = MM(0:N1-N1A-1)
|
||||
end do
|
||||
|
||||
do II = 0, N1A - 1
|
||||
MM(0) = Fyy(II,0)
|
||||
MM(N2A) = Fyy(II,N2A)
|
||||
FF(0:N2A) = F(II,0:N2A)
|
||||
call CreateSpline1 ( CD, CU, N2A + 1, P2, FF, MM, DD, K0, K1, K2 )
|
||||
Fyy(II,0:N2A) = MM(0:N2A)
|
||||
end do
|
||||
|
||||
do II = N1A, N1 - 1
|
||||
MM(0) = Fyy(II,0)
|
||||
MM(N-1) = Fyy(II,N2-1)
|
||||
FF(0:N2-1) = F(II,0:N2-1)
|
||||
call CreateSpline1 ( CD, CU, N2, P2, FF, MM, DD, K0, K1, K2 )
|
||||
Fyy(II,0:N2-1) = MM(0:N2-1)
|
||||
end do
|
||||
|
||||
FF(0:N1-1) = Fyy(0:N1-1,0)
|
||||
call CreateSpline1 ( 3, 3, N1, P1, FF, MM, DD, K0, K1, K2 )
|
||||
Fxxyy(0:N1-1,0) = MM(0:N1-1)
|
||||
|
||||
FF(0:N1A) = Fyy(0:N1A,N2A)
|
||||
call CreateSpline1 ( 3, 3, N1A + 1, P1, FF, MM, DD, K0, K1, K2 )
|
||||
Fxxyy(0:N1A,N2A) = MM(0:N1A)
|
||||
|
||||
FF(0:N1-N1A-1) = Fyy(N1A:N1-1,N2-1 )
|
||||
call CreateSpline1 ( 3, 3, N1-N1A, P1, FF, MM, DD, K0, K1, K2 )
|
||||
Fxxyy(N1A:N1-1,N2-1) = MM(0:N1-N1A-1)
|
||||
|
||||
do II = 1, N1A
|
||||
MM(0) = Fxxyy(II,0)
|
||||
MM(N2A) = Fxxyy(II,N2A)
|
||||
FF(0:N2A) = Fxx(II,0:N2A)
|
||||
call CreateSpline1 ( 2 , 2, N2A + 1, P2, FF, MM, DD, K0, K1, K2 )
|
||||
Fxxyy(II,0:N2A) = MM(0:N2A)
|
||||
end do
|
||||
|
||||
do II = N1A + 1, N1 - 2
|
||||
MM(0) = Fxxyy(II,0)
|
||||
MM(N-1) = Fxxyy(II,N2-1)
|
||||
FF(0:N2-1) = Fxx(II,0:N2-1)
|
||||
call CreateSpline1 ( 2 , 2, N2, P2, FF, MM, DD, K0, K1, K2 )
|
||||
Fxxyy(II,0:N2-1) = MM(0:N2-1)
|
||||
end do
|
||||
|
||||
end subroutine CreateSpline2Ext !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function CalcSpline2_0 ( i, j, X, Y, N1, N2, P1, P2, F, Fxx, Fyy, Fxxyy ) !!!!!!!!!!!
|
||||
integer*4, intent(in) :: i, j, N1, N2
|
||||
real*8, intent(in) :: X, Y
|
||||
real*8, dimension(0:N1-1), intent(in) :: P1
|
||||
real*8, dimension(0:N2-1), intent(in) :: P2
|
||||
real*8, dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy
|
||||
integer*4 :: i1, j1
|
||||
real*8 :: T, Gy_0, Gy_1, Gxxy_0, Gxxy_1
|
||||
!-------------------------------------------------------------------------------------------
|
||||
i1 = i - 1
|
||||
j1 = j - 1
|
||||
T = P2(j) - P2(j1)
|
||||
Gy_0 = ValueSpline1_0 ( Y, P2(j), P2(j1), F(i,j), F(i,j1), Fyy(i,j), Fyy(i,j1), T )
|
||||
Gy_1 = ValueSpline1_0 ( Y, P2(j), P2(j1), F(i1,j), F(i1,j1), Fyy(i1,j), Fyy(i1,j1), T )
|
||||
Gxxy_0 = ValueSpline1_0 ( Y, P2(j), P2(j1), Fxx(i,j), Fxx(i,j1), Fxxyy(i,j), Fxxyy(i,j1), T )
|
||||
Gxxy_1 = ValueSpline1_0 ( Y, P2(j), P2(j1), Fxx(i1,j), Fxx(i1,j1), Fxxyy(i1,j), Fxxyy(i1,j1), T )
|
||||
CalcSpline2_0 = ValueSpline1_0 ( X, P1(i), P1(i1), Gy_0, Gy_1,Gxxy_0, Gxxy_1, P1(i) - P1(i1) )
|
||||
end function CalcSpline2_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine CalcSpline2_1 ( S, Sx1, Sy1, i, j, X, Y, N1, N2, P1, P2, F, Fxx, Fyy, Fxxyy ) !!!
|
||||
real*8, intent(out) :: S, Sx1, Sy1
|
||||
integer*4, intent(in) :: i, j, N1, N2
|
||||
real*8, intent(in) :: X, Y
|
||||
real*8, dimension(0:N1-1), intent(in) :: P1
|
||||
real*8, dimension(0:N2-1), intent(in) :: P2
|
||||
real*8, dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy
|
||||
integer*4 :: i1, j1
|
||||
real*8 :: T, Gy_0, Gy_1, Gxxy_0, Gxxy_1
|
||||
real*8 :: Gyy_0, Gyy_1, Gxxyy_0, Gxxyy_1
|
||||
!-------------------------------------------------------------------------------------------
|
||||
i1 = i - 1
|
||||
j1 = j - 1
|
||||
T = P2(j) - P2(j1)
|
||||
call ValueSpline1_1 ( Gy_0, Gyy_0, Y, P2(j), P2(j1), F(i,j), F(i,j1), Fyy(i,j), Fyy(i,j1), T )
|
||||
call ValueSpline1_1 ( Gy_1, Gyy_1, Y, P2(j), P2(j1), F(i1,j), F(i1,j1), Fyy(i1,j), Fyy(i1,j1), T )
|
||||
call ValueSpline1_1 ( Gxxy_0, Gxxyy_0, Y, P2(j), P2(j1), Fxx(i,j), Fxx(i,j1), Fxxyy(i,j), Fxxyy(i,j1), T )
|
||||
call ValueSpline1_1 ( Gxxy_1, Gxxyy_1, Y, P2(j), P2(j1), Fxx(i1,j), Fxx(i1,j1), Fxxyy(i1,j), Fxxyy(i1,j1), T )
|
||||
call ValueSpline1_1 ( S, Sx1, X, P1(i), P1(i1), Gy_0, Gy_1,Gxxy_0, Gxxy_1, P1(i) - P1(i1) )
|
||||
Sy1 = ValueSpline1_0 ( X, P1(i), P1(i1), Gyy_0, Gyy_1,Gxxyy_0, Gxxyy_1, P1(i) - P1(i1) )
|
||||
end subroutine CalcSpline2_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module Spline2 !********************************************************************************
|
||||
@ -1,62 +0,0 @@
|
||||
program TMDPotGen !*********************************************************************************
|
||||
!
|
||||
! Stand-alone generator of files containing tubular potential data for single-walled CNTs.
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 13.00, 2020
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
use TubePotMono
|
||||
|
||||
implicit none
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Global variables
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4 :: ChiIndM = 10 ! Chirality index m of nanotubes
|
||||
integer*4 :: ChiIndN = 10 ! Chirality index n of nanotubes
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Body
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
TPMStartMode = 0
|
||||
|
||||
! Reading and printing of governing parameters
|
||||
call LoadGoverningParameters ()
|
||||
call PrintGoverningParameters ()
|
||||
|
||||
call TPBInit ()
|
||||
call TPMInit ( ChiIndM, ChiIndN )
|
||||
|
||||
contains !------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine LoadGoverningParameters () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function reads governing parameters from xdt file
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: Fuid, i
|
||||
character*512 :: Msg
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Fuid = OpenFile ( 'TMDPotGen.xdt', 'rt', '' )
|
||||
read ( unit = Fuid, fmt = '(i22)' ) ChiIndM
|
||||
read ( unit = Fuid, fmt = '(i22)' ) ChiIndN
|
||||
call CloseFile ( Fuid )
|
||||
end subroutine LoadGoverningParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine PrintGoverningParameters () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function prints governing parameters to xlg file
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: Fuid, i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Fuid = OpenFile ( 'TMDPotGen.xlg', 'wt', '' )
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) ChiIndM, ' : ChiIndM'
|
||||
write ( unit = Fuid, fmt = '(i22,a)' ) ChiIndN, ' : ChiIndN'
|
||||
call CloseFile ( Fuid )
|
||||
end subroutine PrintGoverningParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end program TMDPotGen !*****************************************************************************
|
||||
@ -1,2 +0,0 @@
|
||||
10 : ChiIndM
|
||||
10 : ChiIndN
|
||||
@ -1,144 +0,0 @@
|
||||
module TPMGeom !************************************************************************************
|
||||
!
|
||||
! Geometry functions for TPM force field.
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 13.00, 2020
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
use TPMLib
|
||||
|
||||
implicit none
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Constants
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4, parameter :: MD_LINES_NONPAR = 0
|
||||
integer*4, parameter :: MD_LINES_PAR = 1
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Global variables
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
! Coordinates of the whole domain
|
||||
real*8 :: DomXmin, DomXmax, DomYmin, DomYmax, DomZmin, DomZmax
|
||||
real*8 :: DomLX, DomLY, DomLZ
|
||||
real*8 :: DomLXhalf, DomLYhalf, DomLZhalf
|
||||
|
||||
! Boundary conditions
|
||||
integer*4 :: BC_X = 0
|
||||
integer*4 :: BC_Y = 0
|
||||
integer*4 :: BC_Z = 0
|
||||
|
||||
! Skin parameter in NBL and related algorithms
|
||||
real*8 :: Rskin = 1.0d+00
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
subroutine ApplyPeriodicBC ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This subroutine changes coordinates of the point according to periodic boundary conditions
|
||||
! it order to make sure that the point is inside the computational cell
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2), intent(inout) :: R
|
||||
!-------------------------------------------------------------------------------------------
|
||||
! These commented lines implement the more general, but less efficient algorithm
|
||||
!if ( BC_X == 1 ) R(0) = R(0) - DomLX * roundint ( R(0) / DomLX )
|
||||
!if ( BC_Y == 1 ) R(1) = R(1) - DomLY * roundint ( R(1) / DomLY )
|
||||
!if ( BC_Z == 1 ) R(2) = R(2) - DomLZ * roundint ( R(2) / DomLZ )
|
||||
if ( BC_X == 1 ) then
|
||||
if ( R(0) .GT. DomLXHalf ) then
|
||||
R(0) = R(0) - DomLX
|
||||
else if ( R(0) .LT. - DomLXHalf ) then
|
||||
R(0) = R(0) + DomLX
|
||||
end if
|
||||
end if
|
||||
if ( BC_Y == 1 ) then
|
||||
if ( R(1) .GT. DomLYHalf ) then
|
||||
R(1) = R(1) - DomLY
|
||||
else if ( R(1) .LT. - DomLYHalf ) then
|
||||
R(1) = R(1) + DomLY
|
||||
end if
|
||||
end if
|
||||
if ( BC_Z == 1 ) then
|
||||
if ( R(2) .GT. DomLZHalf ) then
|
||||
R(2) = R(2) - DomLZ
|
||||
else if ( R(2) .LT. - DomLZHalf ) then
|
||||
R(2) = R(2) + DomLZ
|
||||
end if
|
||||
end if
|
||||
end subroutine ApplyPeriodicBC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine LinePoint ( Displacement, Q, R1, L1, R0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function calculates the point Q of projection of point R0 on line (R1,L1)
|
||||
! Q = R1 + Disaplacement * L1
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(inout) :: Displacement
|
||||
real*8, dimension(0:2), intent(inout) :: Q
|
||||
real*8, dimension(0:2), intent(in) :: R1, L1, R0
|
||||
!--------------------------------------------------------------------------------------------
|
||||
Q = R0 - R1
|
||||
! Here we take into account periodic boundaries
|
||||
call ApplyPeriodicBC ( Q )
|
||||
Displacement = S_V3xV3 ( Q, L1 )
|
||||
Q = R1 + Displacement * L1
|
||||
end subroutine LinePoint !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
integer*4 function LineLine ( H, cosA, D1, D2, L12, R1, L1, R2, L2, Prec ) !!!!!!!!!!!!!!!!!
|
||||
! This function determines the nearest distance H between two lines (R1,L1) and (R2,L2)
|
||||
!-------------------------------------------------------------------------------------------
|
||||
! Input values:
|
||||
! R1, L1, point and direction of line 1
|
||||
! R2, L2, point and direction of line 2
|
||||
! Prec, precision for the case L1 * L2 = 0 (parallel lines)
|
||||
! Return values:
|
||||
! H, minimal distance between lines
|
||||
! cosA, cosine of angle between lines
|
||||
! D1, D2, displacements
|
||||
! L12, unit vector directed along the closes distance
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(inout) :: H, cosA, D1, D2
|
||||
real*8, dimension(0:2), intent(out) :: L12
|
||||
real*8, dimension(0:2), intent(in) :: R1, L1, R2, L2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(in) :: Prec
|
||||
real*8, dimension(0:2) :: Q1, Q2, R
|
||||
real*8 :: C, DD1, DD2, C1, C2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
cosA = S_V3xV3 ( L1, L2 )
|
||||
C = 1.0 - sqr ( cosA )
|
||||
if ( C < Prec ) then ! Lines are parallel to each other
|
||||
LineLine = MD_LINES_PAR
|
||||
return
|
||||
end if
|
||||
LineLine = MD_LINES_NONPAR
|
||||
R = R2 - R1
|
||||
! Here we take into account periodic boundaries
|
||||
call ApplyPeriodicBC ( R )
|
||||
DD1 = S_V3xV3 ( R, L1 )
|
||||
DD2 = S_V3xV3 ( R, L2 )
|
||||
D1 = ( cosA * DD2 - DD1 ) / C
|
||||
D2 = ( DD2 - cosA * DD1 ) / C
|
||||
Q1 = R1 - D1 * L1
|
||||
Q2 = R2 - D2 * L2
|
||||
L12 = Q2 - Q1
|
||||
! Here we take into account periodic boundaries
|
||||
call ApplyPeriodicBC ( L12 )
|
||||
H = S_V3norm3 ( L12 )
|
||||
if ( H < Prec ) then ! Lines intersect each other
|
||||
C1 = signum ( D1 )
|
||||
C2 = signum ( D1 ) * signum ( cosA )
|
||||
Q1 = C1 * L1
|
||||
Q2 = C2 * L2
|
||||
call V3_V3xxV3 ( L12, Q1, Q2 )
|
||||
call V3_ort ( L12 )
|
||||
else ! No intersection
|
||||
L12 = L12 / H
|
||||
end if
|
||||
end function LineLine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module TPMGeom !********************************************************************************
|
||||
@ -1,205 +0,0 @@
|
||||
module TPMLib !*************************************************************************************
|
||||
!
|
||||
! Common constants, types, and functions for TPM force field.
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 13.00, 2020
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
implicit none
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Mathematical constants
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8, parameter :: M_PI_2 = 1.57079632679489661923
|
||||
real*8, parameter :: M_PI = 3.14159265358979323846
|
||||
real*8, parameter :: M_3PI_2 = 4.71238898038468985769
|
||||
real*8, parameter :: M_2PI = 6.28318530717958647692
|
||||
real*8, parameter :: M_PI_180 = 0.017453292519943295769
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Physical unit constants
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8, parameter :: K_AMU = 1.66056E-27 ! a.m.u. (atomic mass unit, Dalton)
|
||||
real*8, parameter :: K_EV = 1.60217646e-19 ! eV (electron-volt)
|
||||
|
||||
real*8, parameter :: K_MDLU = 1.0E-10 ! MD length unit (m)
|
||||
real*8, parameter :: K_MDEU = K_EV ! MD energy unit (J)
|
||||
real*8, parameter :: K_MDMU = K_AMU ! MD mass unit (kg)
|
||||
real*8, parameter :: K_MDFU = K_MDEU / K_MDLU ! MD force unit (N)
|
||||
real*8, parameter :: K_MDCU = K_MDEU / K_MDMU ! MD specific heat unit (J/(kg*K))
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Global variables
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4 :: StdUID = 31
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Simple mathematical functions
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8 function rad ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: X
|
||||
!-------------------------------------------------------------------------------------------
|
||||
rad = X * M_PI_180
|
||||
end function rad !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function sqr ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: X
|
||||
!-------------------------------------------------------------------------------------------
|
||||
sqr = X * X
|
||||
end function sqr !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
integer*4 function signum ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: X
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( X > 0 ) then
|
||||
signum = 1
|
||||
else if ( X < 0 ) then
|
||||
signum = -1
|
||||
else
|
||||
signum = 0
|
||||
end if
|
||||
end function signum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Vector & matrix functions
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8 function S_V3xx ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, dimension(0:2), intent(in) :: V
|
||||
!-------------------------------------------------------------------------------------------
|
||||
S_V3xx = V(0) * V(0) + V(1) * V(1) + V(2) * V(2)
|
||||
end function S_V3xx !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function S_V3xV3 ( V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, dimension(0:2), intent(in) :: V1, V2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
S_V3xV3 = V1(0) * V2(0) + V1(1) * V2(1) + V1(2) * V2(2)
|
||||
end function S_V3xV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function S_V3norm3 ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, dimension(0:2), intent(in) :: V
|
||||
!-------------------------------------------------------------------------------------------
|
||||
S_V3norm3 = dsqrt ( V(0) * V(0) + V(1) * V(1) + V(2) * V(2) )
|
||||
end function S_V3norm3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine V3_ort ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! Vector production
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2), intent(inout) :: V
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: Vabs
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Vabs = S_V3norm3 ( V )
|
||||
V(0) = V(0) / Vabs
|
||||
V(1) = V(1) / Vabs
|
||||
V(2) = V(2) / Vabs
|
||||
end subroutine V3_ort !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine V3_V3xxV3 ( V, V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! Vector production
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2), intent(out) :: V
|
||||
real*8, dimension(0:2), intent(in) :: V1, V2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
V(0) = V1(1) * V2(2) - V1(2) * V2(1)
|
||||
V(1) = V1(2) * V2(0) - V1(0) * V2(2)
|
||||
V(2) = V1(0) * V2(1) - V1(1) * V2(0)
|
||||
end subroutine V3_V3xxV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Handling of spherical and Euler angles
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine RotationMatrix3 ( M, Psi, Tet, Phi ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! Ksi, Tet and Phi are Euler angles
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2,0:2), intent(out) :: M
|
||||
real*8, intent(in) :: Psi, Tet, Phi
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: cK, sK, cT, sT, cP, sP
|
||||
!-------------------------------------------------------------------------------------------
|
||||
cK = dcos ( Psi )
|
||||
sK = dsin ( Psi )
|
||||
cT = dcos ( Tet )
|
||||
sT = dsin ( Tet )
|
||||
cP = dcos ( Phi )
|
||||
sP = dsin ( Phi )
|
||||
M(0,0) = cP * cK - sK * sP * cT
|
||||
M(0,1) = cP * sK + sP * cT * cK
|
||||
M(0,2) = sP * sT
|
||||
M(1,0) = - sP * cK - cP * cT * sK
|
||||
M(1,1) = - sP * sK + cP * cT * cK
|
||||
M(1,2) = cP * sT
|
||||
M(2,0) = sT * sK
|
||||
M(2,1) = - sT * cK
|
||||
M(2,2) = cT
|
||||
end subroutine RotationMatrix3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine EulerAngles ( Psi, Tet, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(out) :: Tet, Psi
|
||||
real*8, dimension(0:2), intent(in) :: L
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Tet = acos ( L(2) )
|
||||
Psi = atan2 ( L(1), L(0) )
|
||||
if ( Psi > M_3PI_2 ) then
|
||||
Psi = Psi - M_3PI_2
|
||||
else
|
||||
Psi = Psi + M_PI_2
|
||||
end if
|
||||
end subroutine EulerAngles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! File input and output
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4 function OpenFile ( Name, Params, Path ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
character*(*), intent(in) :: Name, Params, Path
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: Fuid
|
||||
character*512 :: FullName, Msg, Name1, Action1, Status1, Form1, Position1
|
||||
!-------------------------------------------------------------------------------------------
|
||||
OpenFile = StdUID + 1
|
||||
if ( Params(1:1) == 'r' ) then
|
||||
Action1 = 'read'
|
||||
Status1 = 'old'
|
||||
Position1 = 'rewind'
|
||||
else if ( Params(1:1) == 'w' ) then
|
||||
Action1 = 'write'
|
||||
Status1 = 'replace'
|
||||
Position1 = 'rewind'
|
||||
else if ( Params(1:1) == 'a' ) then
|
||||
Action1 = 'write'
|
||||
Status1 = 'old'
|
||||
Position1 = 'append'
|
||||
endif
|
||||
if ( Params(2:2) == 'b' ) then
|
||||
Form1 = 'binary'
|
||||
else
|
||||
Form1 = 'formatted'
|
||||
endif
|
||||
open ( unit = OpenFile, file = Name, form = Form1, action = Action1, status = Status1, position = Position1 )
|
||||
StdUID = StdUID + 1
|
||||
return
|
||||
end function OpenFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine CloseFile ( Fuid ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
integer*4, intent(inout) :: Fuid
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( Fuid < 0 ) return
|
||||
close ( unit = Fuid )
|
||||
Fuid = -1
|
||||
end subroutine CloseFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module TPMLib !*********************************************************************************
|
||||
|
||||
@ -1,289 +0,0 @@
|
||||
module TubePotBase !********************************************************************************
|
||||
!
|
||||
! Non-Bonded pair interaction potential and transfer functions for atoms composing nanotubes.
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 13.00, 2020
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! This module contains basic parameters for all modules involved into calculations of tubular
|
||||
! potentials.
|
||||
!
|
||||
! It includes definitions of
|
||||
! -- TPBU, Lennard-Jones (12-6) potential;
|
||||
! -- TPBQ, Transfer function,
|
||||
!
|
||||
! All default values are adjusted for non-bonded carbon-carbon interaction in carbon nanotubes.
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
use TPMLib
|
||||
|
||||
implicit none
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Constants
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
! Types of the potential with respect to the breathing mode
|
||||
integer*4, parameter :: TP_POT_MONO_R = 0
|
||||
integer*4, parameter :: TP_POT_POLY_R = 1
|
||||
|
||||
! Maximum number of elements in corresponding tables
|
||||
integer*4, parameter :: TPBNMAX = 2001
|
||||
|
||||
! Numerical constants
|
||||
real*8, parameter :: TPbConstD = 5.196152422706632d+00 ! = 3.0**1.5
|
||||
|
||||
! Mass of C atom
|
||||
real*8, parameter :: TPBMc = 12.0107d+00 ! (Da)
|
||||
|
||||
! Parameters of the Van der Waals interaction between carbon atoms in graphene sheets, see
|
||||
! Stuart S.J., Tutein A.B., Harrison J.A., J. Chem. Phys. 112(14), 2000
|
||||
real*8, parameter :: TPBEcc = 0.00284d+00 ! (eV)
|
||||
real*8, parameter :: TPBScc = 3.4d+00 ! (A)
|
||||
|
||||
! Lattice parameter and surface number density of atoms for a graphene sheet, see
|
||||
! Dresselhaus et al, Carbon 33(7), 1995
|
||||
real*8, parameter :: TPBAcc = 1.421d+00 ! (A)
|
||||
real*8, parameter :: TPBDcc = 4.0d+00 / ( TPBConstD * TPBAcc * TPBAcc ) ! (1/A^2)
|
||||
|
||||
! Specific heat of carbon nanotubes
|
||||
real*8, parameter :: TPBSHcc = 600.0d+00 / K_MDCU ! (eV/(Da*K))
|
||||
|
||||
! Cutoff distances for interactomic potential and transfer function.
|
||||
! Changes in these parameters can result in necessity to change some numerical parameters too.
|
||||
real*8, parameter :: TPBRmincc = 0.001d+00 * TPBScc ! (A)
|
||||
real*8, parameter :: TPBRcutoffcc = 3.0d+00 * TPBScc ! (A)
|
||||
real*8, parameter :: TPBRcutoff1cc = 2.16d+00 * TPBScc ! (A)
|
||||
|
||||
! Parameters of the transfer function for non-bonded interaction between carbon atoms
|
||||
real*8, parameter :: TPBQScc = 7.0d+00 ! (A)
|
||||
real*8, parameter :: TPBQRcutoff1cc = 8.0d+00 ! (A)
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Global variables
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
logical :: TPErrCheck = .true. ! Set to .true. to generate diagnostic and warning messages
|
||||
character*512 :: TPErrMsg = ''
|
||||
|
||||
real*8 :: TPGeomPrec = 1.0d-06 ! Geometric precision, see TPInt
|
||||
integer*4 :: TPPotType = TP_POT_MONO_R ! Type of the potential with respect to the breathing mode
|
||||
|
||||
! Parameters of the interatomic potential and atoms distribution at the nanotube surface
|
||||
|
||||
real*8 :: TPBM = TPBMc ! Mass of an atom (Da)
|
||||
real*8 :: TPBE = TPBEcc ! Depth of the energy well in (12-6) LJ interatomic potential (eV)
|
||||
real*8 :: TPBS = TPBScc ! Sigma parameter of (12-6) LJ interatomic potential (A)
|
||||
real*8 :: TPBD = TPBDcc ! Numerical density of atoms at the tube surface (1/A^2)
|
||||
real*8 :: TPBSH = TPBSHcc ! Specific heat (eV/(Da*K))
|
||||
|
||||
real*8 :: TPBRmin = TPBRmincc ! (A)
|
||||
real*8 :: TPBRcutoff = TPBRcutoffcc ! (A)
|
||||
real*8 :: TPBRcutoff1 = TPBRcutoff1cc ! (A)
|
||||
|
||||
! Physical parameters of the transfer function
|
||||
|
||||
real*8 :: TPBQS = TPBQScc ! Sigma parameter of the transfer function (A)
|
||||
real*8 :: TPBQRcutoff1 = TPBQRcutoff1cc ! (A)
|
||||
|
||||
! Auxiliary variables
|
||||
|
||||
real*8 :: TPBE4, TPBE24, TPBDRcutoff, TPBQDRcutoff
|
||||
real*8 :: TPBQR0 ! Constant-value distance for the transfer function (A)
|
||||
|
||||
! Table of inter-particle potential, force, and transfer function
|
||||
|
||||
integer*4 :: TPBN = TPBNMAX
|
||||
real*8 :: TPBDR
|
||||
real*8, dimension(0:TPBNMAX-1) :: TPBQ
|
||||
real*8, dimension(0:TPBNMAX-1) :: TPBU, TPBdUdR
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
integer*4 function TPBsizeof () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
TPBsizeof = 8 * ( size ( TPBQ ) + size ( TPBU ) + size ( TPBdUdR ) )
|
||||
end function TPBsizeof !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Interpolation
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8 function TPBQInt0 ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: R
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: Z, RR
|
||||
integer*4 :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( R < TPBRmin ) then
|
||||
!call PrintStdLogMsg ( TPErrMsg )
|
||||
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) ': R < Rmin: R=', R, ', Rmin=', TPBRmin
|
||||
!call Error ( 'TPBQInt0', TPErrMsg )
|
||||
elseif ( R > TPBRcutoff ) then
|
||||
TPBQInt0 = 0.0d+00
|
||||
return
|
||||
endif
|
||||
RR = ( R - TPBRmin ) / TPBDR
|
||||
i = int ( RR )
|
||||
RR = RR - i
|
||||
Z = 1.0d+00 - RR
|
||||
TPBQInt0 = TPBQ(i) * Z + TPBQ(i+1) * RR
|
||||
end function TPBQInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function TPBUInt0 ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: R
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: Z, RR
|
||||
integer*4 :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( R < TPBRmin ) then
|
||||
!call PrintStdLogMsg ( TPErrMsg )
|
||||
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) ': R < Rmin: R=', R, ', Rmin=', TPBRmin
|
||||
!call Error ( 'TPBUInt0', TPErrMsg )
|
||||
elseif ( R > TPBRcutoff ) then
|
||||
TPBUInt0 = 0.0d+00
|
||||
return
|
||||
endif
|
||||
RR = ( R - TPBRmin ) / TPBDR
|
||||
i = int ( RR )
|
||||
RR = RR - i
|
||||
Z = 1.0d+00 - RR
|
||||
TPBUInt0 = TPBU(i) * Z + TPBU(i+1) * RR
|
||||
end function TPBUInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine TPBUInt1 ( U, dUdR, R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(out) :: U, dUdR
|
||||
real*8, intent(in) :: R
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: Z, RR
|
||||
integer*4 :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( R < TPBRmin ) then
|
||||
!call PrintStdLogMsg ( TPErrMsg )
|
||||
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) ': R < Rmin: R=', R, ', Rmin=', TPBRmin
|
||||
!call Error ( 'TPBUInt1', TPErrMsg )
|
||||
elseif ( R > TPBRcutoff ) then
|
||||
TPBU = 0.0d+00
|
||||
TPBdUdR = 0.0d+00
|
||||
return
|
||||
endif
|
||||
RR = ( R - TPBRmin ) / TPBDR
|
||||
i = int ( RR )
|
||||
RR = RR - i
|
||||
Z = 1.0d+00 - RR
|
||||
U = TPBU(i) * Z + TPBU(i+1) * RR
|
||||
dUdR = TPBdUdR(i) * Z + TPBdUdR(i+1) * RR
|
||||
end subroutine TPBUInt1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Calculation
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
real*8 function TPBQCalc0 ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: R
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: Z, t, S
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( R > TPBRcutoff ) then
|
||||
TPBQCalc0 = 0.0d+00
|
||||
else if ( R < TPBQR0 ) then
|
||||
TPBQCalc0 = 1.0d+00
|
||||
else
|
||||
Z = TPBQS / R
|
||||
Z = Z * Z * Z
|
||||
Z = Z * Z
|
||||
TPBQCalc0 = 4.0d+00 * ( 1.0d+00 - Z ) * Z
|
||||
if ( R > TPBQRcutoff1 ) then
|
||||
t = ( R - TPBQRcutoff1 ) / TPBQDRcutoff
|
||||
S = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
|
||||
TPBQCalc0 = TPBQCalc0 * S
|
||||
endif
|
||||
endif
|
||||
end function TPBQCalc0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
real*8 function TPBUCalc0 ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: R
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: Z, t, S
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( R > TPBRcutoff ) then
|
||||
TPBUCalc0 = 0.0d+00
|
||||
else
|
||||
Z = TPBS / R
|
||||
Z = Z * Z * Z
|
||||
Z = Z * Z
|
||||
TPBUCalc0 = TPBE4 * ( Z - 1.0d+00 ) * Z
|
||||
if ( R > TPBRcutoff1 ) then
|
||||
t = ( R - TPBRcutoff1 ) / TPBDRcutoff
|
||||
S = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
|
||||
TPBUCalc0 = TPBUCalc0 * S
|
||||
endif
|
||||
endif
|
||||
end function TPBUCalc0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine TPBUCalc1 ( U, dUdR, R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(out) :: U, dUdR
|
||||
real*8, intent(in) :: R
|
||||
real*8 :: Z, t, S, dSdR
|
||||
!-------------------------------------------------------------------------------------------
|
||||
if ( R > TPBRcutoff ) then
|
||||
U = 0.0d+00
|
||||
dUdR = 0.0d+00
|
||||
else
|
||||
Z = TPBS / R
|
||||
Z = Z * Z * Z
|
||||
Z = Z * Z
|
||||
U = TPBE4 * ( Z - 1.0d+00 ) * Z
|
||||
dUdR = TPBE24 * ( 2.0d+00 * Z - 1.0d+00 ) * Z / R
|
||||
if ( R > TPBRcutoff1 ) then
|
||||
t = ( R - TPBRcutoff1 ) / TPBDRcutoff
|
||||
S = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
|
||||
dSdR = 6.0d+00 * t * ( t - 1.0d+00 ) / TPBDRcutoff
|
||||
dUdR = dUdR * S + U * dSdR
|
||||
U = U * S
|
||||
endif
|
||||
endif
|
||||
end subroutine TPBUCalc1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine TPBSegmentForces ( F1, F2, F, M, Laxis, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, dimension(0:2), intent(out) :: F1, F2
|
||||
real*8, dimension(0:2), intent(in) :: F, M, Laxis
|
||||
real*8, intent(in) :: L
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2) :: FF, MM, FFF
|
||||
!-------------------------------------------------------------------------------------------
|
||||
FF = 0.5d+00 * F
|
||||
MM = M / L
|
||||
call V3_V3xxV3 ( FFF, MM, Laxis )
|
||||
F1 = FF - FFF
|
||||
F2 = FF + FFF
|
||||
end subroutine TPBSegmentForces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Initialization
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine TPBInit () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8 :: R
|
||||
integer*4 :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
TPBE4 = 4.0d+00 * TPBE
|
||||
TPBE24 = - 24.0d+00 * TPBE
|
||||
TPBDRcutoff = TPBRcutoff - TPBRcutoff1
|
||||
TPBQDRcutoff = TPBRcutoff - TPBQRcutoff1
|
||||
TPBQR0 = TPBQS * 2.0d+00 ** ( 1.0d+00 / 6.0d+00 )
|
||||
TPBDR = ( TPBRcutoff - TPBRmin ) / ( TPBN - 1 )
|
||||
R = TPBRmin
|
||||
do i = 0, TPBN - 1
|
||||
TPBQ(i) = TPBQCalc0 ( R )
|
||||
call TPBUCalc1 ( TPBU(i), TPBdUdR(i), R )
|
||||
R = R + TPBDR
|
||||
enddo
|
||||
end subroutine TPBInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module TubePotBase !****************************************************************************
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,428 +0,0 @@
|
||||
module TubePotTrue !********************************************************************************
|
||||
!
|
||||
! True tubular potential and transfer function
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! This module implements calculation of true potential and transfer functions for interaction
|
||||
! between two cylinder segments of nanotubes by direct integration over the surfaces of both
|
||||
! segments.
|
||||
!
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
!
|
||||
! Intel Fortran
|
||||
!
|
||||
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 13.00, 2020
|
||||
!
|
||||
!***************************************************************************************************
|
||||
|
||||
use TPMGeom
|
||||
use TubePotBase
|
||||
|
||||
implicit none
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Constants
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4, parameter :: TPTNXMAX = 257
|
||||
integer*4, parameter :: TPTNEMAX = 128
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Types
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
type TPTSEG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8 :: X, Y, Z
|
||||
real*8 :: Psi, Theta, Phi ! Euler's angles
|
||||
real*8 :: R ! Segment radius
|
||||
real*8 :: L ! Segment length
|
||||
integer*4 :: NX, NE ! Number of nodes for numerical integration
|
||||
real*8 :: DX, DE ! Spacings
|
||||
real*8, dimension(0:2,0:2) :: M ! Transformation matrix
|
||||
real*8, dimension(0:TPTNXMAX-1,0:TPTNXMAX-1,0:2) :: Rtab! Node coordinates
|
||||
end type TPTSEG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Global variables
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
type(TPTSEG) :: TPTSeg1, TPTSeg2 ! Two segments
|
||||
|
||||
contains !******************************************************************************************
|
||||
|
||||
subroutine TPTSegAxisVector ( S, Laxis ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
type(TPTSEG), intent(in) :: S
|
||||
real*8, dimension(0:2), intent(out) :: Laxis
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Laxis(0:2) = S%M(2,0:2)
|
||||
end subroutine TPTSegAxisVector !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine TPTSegRadVector ( S, Lrad, Eps ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
type(TPTSEG), intent(in) :: S
|
||||
real*8, dimension(0:2), intent(out) :: Lrad
|
||||
real*8, intent(in) :: Eps
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: Ce, Se
|
||||
!-------------------------------------------------------------------------------------------
|
||||
Ce = cos ( Eps )
|
||||
Se = sin ( Eps )
|
||||
Lrad(0) = Ce * S%M(0,0) + Se * S%M(1,0)
|
||||
Lrad(1) = Ce * S%M(0,1) + Se * S%M(1,1)
|
||||
Lrad(2) = Ce * S%M(0,2) + Se * S%M(1,2)
|
||||
end subroutine TPTSegRadVector !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine TPTRadiusVector ( S, R, X, Eps ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
type(TPTSEG), intent(in) :: S
|
||||
real*8, dimension(0:2), intent(out) :: R
|
||||
real*8, intent(in) :: X, Eps
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2) :: Laxis, Lrad
|
||||
!-------------------------------------------------------------------------------------------
|
||||
call TPTSegAxisVector ( S, Laxis )
|
||||
call TPTSegRadVector ( S, Lrad, Eps )
|
||||
R(0) = S%X + X * Laxis(0) + S%R * Lrad(0)
|
||||
R(1) = S%Y + X * Laxis(1) + S%R * Lrad(1)
|
||||
R(2) = S%Z + X * Laxis(2) + S%R * Lrad(2)
|
||||
end subroutine TPTRadiusVector !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine TPTCalcSegNodeTable ( S ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
type(TPTSEG), intent(inout) :: S
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: X, Eps
|
||||
integer*4 :: i, j
|
||||
!-------------------------------------------------------------------------------------------
|
||||
X = - S%L / 2.0
|
||||
call RotationMatrix3 ( S%M, S%Psi, S%Theta, S%Phi )
|
||||
do i = 0, S%NX - 1
|
||||
Eps = 0.0d+00
|
||||
do j = 0, S%NE - 1
|
||||
call TPTRadiusVector ( S, S%Rtab(i,j,0:2), X, Eps )
|
||||
Eps = Eps + S%DE
|
||||
end do
|
||||
X = X + S%DX
|
||||
end do
|
||||
end subroutine TPTCalcSegNodeTable !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine TPTSetSegPosition1 ( S, Rcenter, Laxis, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
type(TPTSEG), intent(inout) :: S
|
||||
real*8, dimension(0:2), intent(in) :: Rcenter, Laxis
|
||||
real*8, intent(in) :: L
|
||||
!-------------------------------------------------------------------------------------------
|
||||
S%L = L
|
||||
S%DX = L / ( S%NX - 1 )
|
||||
call EulerAngles ( S%Psi, S%Theta, Laxis )
|
||||
S%Phi= 0.0d+00
|
||||
S%X = Rcenter(0)
|
||||
S%Y = Rcenter(1)
|
||||
S%Z = Rcenter(2)
|
||||
call TPTCalcSegNodeTable ( S )
|
||||
end subroutine TPTSetSegPosition1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine TPTSetSegPosition2 ( S, R1, R2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
type(TPTSEG), intent(inout) :: S
|
||||
real*8, dimension(0:2), intent(in) :: R1, R2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2) :: R, Laxis
|
||||
real*8 :: L
|
||||
!-------------------------------------------------------------------------------------------
|
||||
R = 0.5 * ( R1 + R2 )
|
||||
Laxis = R2 - R1
|
||||
L = S_V3norm3 ( Laxis )
|
||||
Laxis = Laxis / L
|
||||
call TPTSetSegPosition1 ( S, R, Laxis, L )
|
||||
end subroutine TPTSetSegPosition2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
integer*4 function TPTCheckIntersection ( S1, S2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
type(TPTSEG), intent(in) :: S1, S2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: i, j
|
||||
real*8 :: L1, L2, Displacement, D
|
||||
real*8, dimension(0:2) :: Laxis, Q, R
|
||||
!-------------------------------------------------------------------------------------------
|
||||
L2 = S1%L / 2.0
|
||||
L1 = - L2
|
||||
call TPTSegAxisVector ( S1, Laxis )
|
||||
R(0) = S1%X
|
||||
R(1) = S1%Y
|
||||
R(2) = S1%Z
|
||||
do i = 0, S2%NX - 1
|
||||
do j = 0, S2%NE - 1
|
||||
call LinePoint ( Displacement, Q, R, Laxis, S2%Rtab(i,j,0:2) )
|
||||
D = sqrt ( sqr ( Q(0) - S2%Rtab(i,j,0) ) + sqr ( Q(1) - S2%Rtab(i,j,1) ) &
|
||||
+ sqr ( Q(2) - S2%Rtab(i,j,2) ) )
|
||||
if ( Displacement > L1 .and. Displacement < L2 .and. D < S1%R ) then
|
||||
TPTCheckIntersection = 1
|
||||
return
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
TPTCheckIntersection = 0
|
||||
end function TPTCheckIntersection !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
integer*4 function TPTCalcPointRange ( S, Xmin, Xmax, Re ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
type(TPTSEG), intent(in) :: S
|
||||
real*8, intent(out) :: Xmin, Xmax
|
||||
real*8, dimension(0:2), intent(in) :: Re
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8 :: Displacement, Distance
|
||||
real*8, dimension(0:2) :: Laxis, Q, R
|
||||
!-------------------------------------------------------------------------------------------
|
||||
call TPTSegAxisVector ( S, Laxis )
|
||||
R(0) = S%X
|
||||
R(1) = S%Y
|
||||
R(2) = S%Z
|
||||
call LinePoint ( Displacement, Q, R, Laxis, Re )
|
||||
Distance = sqrt ( sqr ( Q(0) - Re(0) ) + sqr ( Q(1) - Re(1) ) + sqr ( Q(2) - Re(2) ) ) - S%R
|
||||
if ( TPBRcutoff < Distance ) then
|
||||
Xmin = 0.0d+00
|
||||
Xmax = 0.0d+00
|
||||
TPTCalcPointRange = 0
|
||||
return
|
||||
end if
|
||||
Distance = sqrt ( TPBRcutoff * TPBRcutoff - Distance * Distance )
|
||||
Xmin = Displacement - Distance
|
||||
Xmax = Displacement + Distance
|
||||
TPTCalcPointRange = 1
|
||||
end function TPTCalcPointRange !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
subroutine TPTGetEnds ( R1_1, R1_2, R2_1, R2_2, X1_1, X1_2, X2_1, X2_2, H, A ) !!!!!!!!!!!!!
|
||||
real*8, dimension(0:2), intent(out) :: R1_1, R1_2, R2_1, R2_2
|
||||
real*8, intent(in) :: X1_1, X1_2, X2_1, X2_2, H, A
|
||||
!-------------------------------------------------------------------------------------------
|
||||
R1_1(0) = 0.0d+00
|
||||
R1_1(1) = 0.0d+00
|
||||
R1_1(2) = X1_1
|
||||
R1_2(0) = 0.0d+00
|
||||
R1_2(1) = 0.0d+00
|
||||
R1_2(2) = X1_2
|
||||
R2_1(0) = H
|
||||
R2_1(1) = - X2_1 * sin ( A )
|
||||
R2_1(2) = X2_1 * cos ( A )
|
||||
R2_2(0) = H
|
||||
R2_2(1) = - X2_2 * sin ( A )
|
||||
R2_2(2) = X2_2 * cos ( A )
|
||||
end subroutine TPTGetEnds !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Tubular potential
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
integer*4 function TPTPointPotential ( Q, U, F, R, S ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function returns the potential U and force F applied to the atom in position R and
|
||||
! produced by the segment S.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(out) :: Q, U
|
||||
real*8, dimension(0:2), intent(out) :: F
|
||||
real*8, dimension(0:2), intent(in) :: R
|
||||
type(TPTSEG), intent(in) :: S
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: i, j
|
||||
real*8, dimension(0:2) :: RR, FF
|
||||
real*8 :: QQ, UU, UUU, FFF, Rabs
|
||||
real*8 :: Coeff, Xmin, Xmax, X
|
||||
!-------------------------------------------------------------------------------------------
|
||||
TPTPointPotential = 0
|
||||
Q = 0.0d+00
|
||||
U = 0.0d+00
|
||||
F = 0.0d+00
|
||||
if ( TPTCalcPointRange ( S, Xmin, Xmax, R ) == 0 ) return
|
||||
X = - S%L / 2.0
|
||||
do i = 0, S%NX - 1
|
||||
if ( X > Xmin .and. X < Xmax ) then
|
||||
QQ = 0.0d+00
|
||||
UU = 0.0d+00
|
||||
FF = 0.0d+00
|
||||
do j = 0, S%NE - 1
|
||||
RR(0:2) = S%Rtab(i,j,0:2) - R(0:2)
|
||||
Rabs = S_V3norm3 ( RR )
|
||||
if ( Rabs < TPBRcutoff ) then
|
||||
QQ = QQ + TPBQCalc0 ( Rabs )
|
||||
call TPBUCalc1 ( UUU, FFF, Rabs )
|
||||
UU = UU + UUU
|
||||
FFF = FFF / Rabs
|
||||
FF = FF + FFF * RR
|
||||
TPTPointPotential = 1
|
||||
end if
|
||||
end do
|
||||
if ( i == 0 .or. i == S%NX - 1 ) then
|
||||
Q = Q + 0.5d+00 * QQ
|
||||
U = U + 0.5d+00 * UU
|
||||
F = F + 0.5d+00 * FF
|
||||
else
|
||||
Q = Q + QQ
|
||||
U = U + UU
|
||||
F = F + FF
|
||||
end if
|
||||
end if
|
||||
X = X + S%DX
|
||||
end do
|
||||
Coeff = TPBD * S%DX * S%R * S%DE
|
||||
Q = Q * S%DX * S%R * S%DE
|
||||
U = U * Coeff
|
||||
F = F * Coeff
|
||||
end function TPTPointPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
integer*4 function TPTSectionPotential ( Q, U, F, M, S, i, Ssource ) !!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function returns the potential U, force F and torque M produced by the segment Ssource
|
||||
! and applied to the i-th circular cross-section of the segment S.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(out) :: Q, U
|
||||
real*8, dimension(0:2), intent(out) :: F, M
|
||||
type(TPTSEG), intent(in) :: S, Ssource
|
||||
integer*4, intent(in) :: i
|
||||
!-------------------------------------------------------------------------------------------
|
||||
integer*4 :: j
|
||||
real*8, dimension(0:2) :: R, Fp, Mp, Lrad
|
||||
real*8 :: Qp, Up, Eps
|
||||
real*8 :: Coeff
|
||||
!-------------------------------------------------------------------------------------------
|
||||
TPTSectionPotential = 0
|
||||
Q = 0.0d+00
|
||||
U = 0.0d+00
|
||||
F = 0.0d+00
|
||||
M = 0.0d+00
|
||||
Eps = 0.0d+00
|
||||
do j = 0, S%NE - 1
|
||||
call TPTSegRadVector ( S, Lrad, Eps )
|
||||
if ( TPTPointPotential ( Qp, Up, Fp, S%Rtab(i,j,0:2), Ssource ) == 1 ) then
|
||||
Q = Q + Qp
|
||||
U = U + Up
|
||||
F = F + Fp
|
||||
R(0) = S%Rtab(i,j,0) - S%X
|
||||
R(1) = S%Rtab(i,j,1) - S%Y
|
||||
R(2) = S%Rtab(i,j,2) - S%Z
|
||||
call V3_V3xxV3 ( Mp, R, Fp )
|
||||
M = M + Mp
|
||||
TPTSectionPotential = 1
|
||||
end if
|
||||
Eps = Eps + S%DE
|
||||
end do
|
||||
Coeff = TPBD * S%R * S%DE
|
||||
Q = Q * S%R * S%DE
|
||||
U = U * Coeff
|
||||
F = F * Coeff
|
||||
M = M * Coeff
|
||||
end function TPTSectionPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
integer*4 function TPTSegmentPotential ( Q, U, F, M, S, Ssource ) !!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! This function returns the potential U, force F and torque M produced by the segment
|
||||
! Ssource and applied to the segment S.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(out) :: Q, U
|
||||
real*8, dimension(0:2), intent(out) :: F, M
|
||||
type(TPTSEG), intent(in) :: S, Ssource
|
||||
integer*4 :: i
|
||||
real*8, dimension(0:2) :: Fc, Mc
|
||||
real*8 :: Qc, Uc
|
||||
!-------------------------------------------------------------------------------------------
|
||||
TPTSegmentPotential = 0
|
||||
Q = 0.0d+00
|
||||
U = 0.0d+00
|
||||
F = 0.0d+00
|
||||
M = 0.0d+00
|
||||
if ( TPTCheckIntersection ( S, Ssource ) == 1 ) then
|
||||
TPTSegmentPotential = 2
|
||||
return
|
||||
end if
|
||||
do i = 0, S%NX - 1
|
||||
if ( TPTSectionPotential ( Qc, Uc, Fc, Mc, S, i, Ssource ) == 1 ) then
|
||||
if ( i == 0 .or. i == S%NX - 1 ) then
|
||||
Q = Q + 0.5d+00 * Qc
|
||||
U = U + 0.5d+00 * Uc
|
||||
F = F + 0.5d+00 * Fc
|
||||
M = M + 0.5d+00 * Mc
|
||||
else
|
||||
Q = Q + Qc
|
||||
U = U + Uc
|
||||
F = F + Fc
|
||||
M = M + Mc
|
||||
end if
|
||||
TPTSegmentPotential = 1
|
||||
end if
|
||||
end do
|
||||
Q = Q * S%DX
|
||||
U = U * S%DX
|
||||
F = F * S%DX
|
||||
M = M * S%DX
|
||||
end function TPTSegmentPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Forces
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine TPTSegmentForces ( F1, F2, F, M, Laxis, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, dimension(0:2), intent(out) :: F1, F2
|
||||
real*8, dimension(0:2), intent(in) :: F, M, Laxis
|
||||
real*8, intent(in) :: L
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2) :: MM, FF, FFF
|
||||
!-------------------------------------------------------------------------------------------
|
||||
FF = 0.5d+00 * F
|
||||
MM = M / L
|
||||
call V3_V3xxV3 ( FFF, MM, Laxis )
|
||||
F1 = FF - FFF
|
||||
F2 = FF + FFF
|
||||
end subroutine TPTSegmentForces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
integer*4 function TPTInteractionF ( Q, U, F1_1, F1_2, F2_1, F2_2, R1_1, R1_2, R2_1, R2_2 )
|
||||
! This function returns the potential and forces applied to the ends of segments.
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, intent(out) :: Q, U
|
||||
real*8, dimension(0:2), intent(out) :: F1_1, F1_2, F2_1, F2_2
|
||||
real*8, dimension(0:2), intent(in) :: R1_1, R1_2, R2_1, R2_2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
real*8, dimension(0:2) :: R1, R2, Laxis1, Laxis2, DR, F1, M1, F2, M2
|
||||
real*8 :: L1, L2
|
||||
!-------------------------------------------------------------------------------------------
|
||||
R1 = 0.5 * ( R1_1 + R1_2 )
|
||||
R2 = 0.5 * ( R2_1 + R2_2 )
|
||||
Laxis1 = R1_2 - R1_1
|
||||
Laxis2 = R2_2 - R2_1
|
||||
L1 = S_V3norm3 ( Laxis1 )
|
||||
L2 = S_V3norm3 ( Laxis2 )
|
||||
Laxis1 = Laxis1 / L1
|
||||
Laxis2 = Laxis2 / L2
|
||||
DR = R2 - R1
|
||||
call TPTSetSegPosition1 ( TPTSeg1, R1, Laxis1, L1 )
|
||||
call TPTSetSegPosition1 ( TPTSeg2, R2, Laxis2, L2 )
|
||||
TPTInteractionF = TPTSegmentPotential ( Q, U, F1, M1, TPTSeg1, TPTSeg2 )
|
||||
if ( TPTInteractionF .ne. 1 ) return
|
||||
call V3_V3xxV3 ( M2, DR, F1 )
|
||||
F2 = - F1
|
||||
M2 = - M1 - M2
|
||||
call TPTSegmentForces ( F1_1, F1_2, F1, M1, Laxis1, L1 )
|
||||
call TPTSegmentForces ( F2_1, F2_2, F2, M2, Laxis2, L2 )
|
||||
end function TPTInteractionF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
! Initialization
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
|
||||
subroutine TPTInit ( R1, R2, NX, NE ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
real*8, intent(in) :: R1, R2
|
||||
integer*4, intent(in) :: NX, NE
|
||||
!-------------------------------------------------------------------------------------------
|
||||
TPTSeg1%X = 0.0d+00
|
||||
TPTSeg1%Y = 0.0d+00
|
||||
TPTSeg1%Z = 0.0d+00
|
||||
TPTSeg1%Psi = 0.0d+00
|
||||
TPTSeg1%Theta = 0.0d+00
|
||||
TPTSeg1%Phi = 0.0d+00
|
||||
TPTSeg1%R = R1
|
||||
TPTSeg1%NX = NX
|
||||
TPTSeg1%NE = NE
|
||||
TPTSeg1%DE = M_2PI / NE
|
||||
TPTSeg2%X = 0.0d+00
|
||||
TPTSeg2%Y = 0.0d+00
|
||||
TPTSeg2%Z = 0.0d+00
|
||||
TPTSeg2%Psi = 0.0d+00
|
||||
TPTSeg2%Theta = 0.0d+00
|
||||
TPTSeg2%Phi = 0.0d+00
|
||||
TPTSeg2%R = R2
|
||||
TPTSeg2%NX = NX
|
||||
TPTSeg2%NE = NE
|
||||
TPTSeg2%DE = M_2PI / NE
|
||||
end subroutine TPTInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end module TubePotTrue !****************************************************************************
|
||||
@ -1,247 +0,0 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
|
||||
Contributing author: Maxim Shugaev (UVA), mvs9t@virginia.edu
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <iostream>
|
||||
#include <cstdlib>
|
||||
#include <fstream>
|
||||
#include <string>
|
||||
#include <string.h>
|
||||
#include <vector>
|
||||
#include <array>
|
||||
#include <regex>
|
||||
#include <string.h>
|
||||
#include <cmath>
|
||||
//#include <filesystem>
|
||||
|
||||
static const std::string data_file0 = "system.init";
|
||||
static const std::string data_dump0 = "config.dump";
|
||||
static const std::string out_dir0 = "out";
|
||||
|
||||
struct Particle {
|
||||
double x, y, z, vx, vy, vz, Es, Eb, Et, Ep, Ek;
|
||||
char type, nx, ny, nz;
|
||||
};
|
||||
|
||||
class Lamps_base {
|
||||
public:
|
||||
Lamps_base() = default;
|
||||
virtual ~Lamps_base() = default;
|
||||
|
||||
int open(const std::string& filename);
|
||||
int next(); //get next snapshot from the opened file
|
||||
virtual int write(const std::string& filename) const = 0;
|
||||
|
||||
inline double get_X1() const { return X1; };
|
||||
inline double get_X2() const { return X2; };
|
||||
inline double get_Y1() const { return Y1; };
|
||||
inline double get_Y2() const { return Y2; };
|
||||
inline double get_Z1() const { return Z1; };
|
||||
inline double get_Z2() const { return Z2; };
|
||||
inline int get_Natoms() const { return Natoms; };
|
||||
inline int get_Nsteps() const { return Nsteps; };
|
||||
inline int is_open() const { return open_stat; };
|
||||
inline const Particle& get(int i) const { return particles[i]; };
|
||||
inline Particle& get(int i) { return particles[i]; };
|
||||
|
||||
protected:
|
||||
virtual int load() = 0;
|
||||
|
||||
int Nsteps, Natoms, open_stat;
|
||||
double X1, X2, Y1, Y2, Z1, Z2;
|
||||
std::vector<Particle> particles;
|
||||
std::ifstream in;
|
||||
};
|
||||
|
||||
class Lamps_dump : public Lamps_base {
|
||||
public:
|
||||
Lamps_dump() = default;
|
||||
~Lamps_dump() = default;
|
||||
|
||||
virtual int write(const std::string& filename) const override;
|
||||
private:
|
||||
virtual int load() override;
|
||||
};
|
||||
|
||||
int Lamps_base::open(const std::string& filename) {
|
||||
in.open(filename); if (!in.is_open()) return EXIT_FAILURE;
|
||||
return load();
|
||||
}
|
||||
|
||||
int Lamps_base::next() {
|
||||
return load();
|
||||
}
|
||||
|
||||
int Lamps_dump::write(const std::string& filename) const {
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
|
||||
int Lamps_dump::load() {
|
||||
std::string inbuf; char* tmp_cptr;
|
||||
open_stat = 0;
|
||||
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
Nsteps = std::stoi(inbuf);
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
Natoms = std::stoi(inbuf);
|
||||
particles.resize(Natoms);
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
X1 = strtof(inbuf.c_str(), &tmp_cptr);
|
||||
X2 = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
Y1 = strtof(inbuf.c_str(), &tmp_cptr);
|
||||
Y2 = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
Z1 = strtof(inbuf.c_str(), &tmp_cptr);
|
||||
Z2 = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
for (int i = 0; i < Natoms; i++) {
|
||||
if (!getline(in, inbuf)) return EXIT_FAILURE;
|
||||
int id = strtol(inbuf.c_str(), &tmp_cptr, 10) - 1; // modify based on a particular file format
|
||||
particles[id].type = static_cast<char>(strtol(tmp_cptr + 1, &tmp_cptr, 10));
|
||||
particles[id].x = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
particles[id].y = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
particles[id].z = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
particles[id].Es = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
particles[id].Eb = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
particles[id].Et = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
particles[id].Ep = particles[id].Es + particles[id].Eb + particles[id].Et;
|
||||
particles[id].Ek = strtof(tmp_cptr + 1, &tmp_cptr);
|
||||
}
|
||||
open_stat = true;
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char* argv[]) {
|
||||
std::string data_file = (argc > 1) ? argv[1] : data_file0;
|
||||
std::string data_dump = (argc > 2) ? argv[2] : data_dump0;
|
||||
std::string out_dir = (argc > 3) ? argv[3] : out_dir0;
|
||||
//std::filesystem::remove_all(out_dir);
|
||||
//std::filesystem::create_directories(out_dir);
|
||||
|
||||
//list of bonds
|
||||
std::ifstream in(data_file);
|
||||
if (!in.is_open()) {
|
||||
std::cout << "cannot open " << data_file << std::endl;
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
std::string buf;
|
||||
std::string atoms_l = "Atoms";
|
||||
while (std::getline(in, buf)){
|
||||
if (buf == atoms_l) break;
|
||||
if (in.eof()) return EXIT_FAILURE;
|
||||
}
|
||||
std::getline(in, buf);
|
||||
char* tmp_cptr;
|
||||
std::vector<std::array<int, 2>> bonds;
|
||||
while (std::getline(in, buf)) {
|
||||
if (in.eof() || buf.size() == 0) break;
|
||||
int idx = strtol(buf.c_str(), &tmp_cptr, 10);
|
||||
int m_idx = strtol(tmp_cptr + 1, &tmp_cptr, 10);
|
||||
int type = strtol(tmp_cptr + 1, &tmp_cptr, 10);
|
||||
int id1 = strtol(tmp_cptr + 1, &tmp_cptr, 10);
|
||||
int id2 = strtol(tmp_cptr + 1, &tmp_cptr, 10);
|
||||
if(id1 >= 0 && id2 >= 0) bonds.push_back({id1 - 1, id2 - 1});
|
||||
}
|
||||
|
||||
//dump
|
||||
Lamps_dump dump;
|
||||
dump.open(data_dump);
|
||||
if (!dump.is_open()) {
|
||||
std::cout << "cannot open " << data_dump << std::endl;
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
double Lx = dump.get_X2() - dump.get_X1();
|
||||
double Ly = dump.get_Y2() - dump.get_Y1();
|
||||
double Lz = dump.get_Z2() - dump.get_Z1();
|
||||
while (1) {
|
||||
std::ofstream out(out_dir + "/cnt" + std::to_string(dump.get_Nsteps()) + ".vtk");
|
||||
if (!out.is_open()) {
|
||||
std::cout << "cannot create " << out_dir + "/cnt" + std::to_string(dump.get_Nsteps()) + ".vtk" << std::endl;
|
||||
std::cout << "create the output directory \"" << out_dir << "\" manually" << std::endl;
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
out << "# vtk DataFile Version 3.0\n# \nASCII\n\nDATASET UNSTRUCTURED_GRID\n";
|
||||
out << "POINTS " << dump.get_Natoms() << " float\n";
|
||||
for (int i = 0; i < dump.get_Natoms(); i++) {
|
||||
out << dump.get(i).x << " " << dump.get(i).y << " " << dump.get(i).z << " " << "\n";
|
||||
}
|
||||
|
||||
int bond_count = 0;
|
||||
for (int i = 0; i < bonds.size(); i++) {
|
||||
double f1 = std::fabs(dump.get(bonds[i][0]).x - dump.get(bonds[i][1]).x);
|
||||
double f2 = std::fabs(dump.get(bonds[i][0]).y - dump.get(bonds[i][1]).y);
|
||||
double f3 = std::fabs(dump.get(bonds[i][0]).z - dump.get(bonds[i][1]).z);
|
||||
if ((std::fabs(dump.get(bonds[i][0]).x - dump.get(bonds[i][1]).x) < 0.5*Lx)
|
||||
&& (std::fabs(dump.get(bonds[i][0]).y - dump.get(bonds[i][1]).y) < 0.5*Ly)
|
||||
&& (std::fabs(dump.get(bonds[i][0]).z - dump.get(bonds[i][1]).z) < 0.5*Lz))
|
||||
bond_count++;
|
||||
}
|
||||
out << "\nCELLS " << bond_count << " " << 3*bond_count << "\n";
|
||||
for (int i = 0; i < bonds.size(); i++) {
|
||||
if ((std::fabs(dump.get(bonds[i][0]).x - dump.get(bonds[i][1]).x) < 0.5 * Lx)
|
||||
&& (std::fabs(dump.get(bonds[i][0]).y - dump.get(bonds[i][1]).y) < 0.5 * Ly)
|
||||
&& (std::fabs(dump.get(bonds[i][0]).z - dump.get(bonds[i][1]).z) < 0.5 * Lz))
|
||||
out << "2 " << bonds[i][0] << " " << bonds[i][1] << " " << "\n";
|
||||
}
|
||||
|
||||
out << "\nCELL_TYPES " << bond_count << "\n";
|
||||
for (int i = 0; i < bond_count; i++) {
|
||||
out << "4\n";
|
||||
}
|
||||
|
||||
out << "\nPOINT_DATA " << dump.get_Natoms() << "\n";
|
||||
out << "SCALARS Ep float 1\n";
|
||||
out << "LOOKUP_TABLE default\n";
|
||||
for (int i = 0; i < dump.get_Natoms(); i++) {
|
||||
out << dump.get(i).Ep << "\n";
|
||||
}
|
||||
|
||||
out << "\nSCALARS Ek float 1\n";
|
||||
out << "LOOKUP_TABLE default\n";
|
||||
for (int i = 0; i < dump.get_Natoms(); i++) {
|
||||
out << dump.get(i).Ek << "\n";
|
||||
}
|
||||
|
||||
out << "\nSCALARS Es float 1\n";
|
||||
out << "LOOKUP_TABLE default\n";
|
||||
for (int i = 0; i < dump.get_Natoms(); i++) {
|
||||
out << dump.get(i).Es << "\n";
|
||||
}
|
||||
|
||||
out << "\nSCALARS Eb float 1\n";
|
||||
out << "LOOKUP_TABLE default\n";
|
||||
for (int i = 0; i < dump.get_Natoms(); i++) {
|
||||
out << dump.get(i).Eb << "\n";
|
||||
}
|
||||
|
||||
out << "\nSCALARS Et float 1\n";
|
||||
out << "LOOKUP_TABLE default\n";
|
||||
for (int i = 0; i < dump.get_Natoms(); i++) {
|
||||
out << dump.get(i).Et << "\n";
|
||||
}
|
||||
|
||||
if (dump.next() != EXIT_SUCCESS) break;
|
||||
}
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
Reference in New Issue
Block a user