Commit 0023c848 authored by Jony Castagna's avatar Jony Castagna
Browse files

Revert "Merge branch 'dlmeso_tetrahedral' into 'master'"

This reverts merge request !64
parent 5cb0ab5c
......@@ -62,7 +62,6 @@ The following modules connected to the DL_MESO_DPD code have been produced so fa
One species - starting as diamond cubic lattice
0 1
4.0000000000 0.0000000000 0.0000000000
0.0000000000 4.0000000000 0.0000000000
0.0000000000 0.0000000000 4.0000000000
A 1
0.0 0.0 0.0
A 2
0.0 2.0 2.0
A 3
2.0 0.0 2.0
A 4
2.0 2.0 0.0
A 5
3.0 3.0 3.0
A 6
3.0 1.0 1.0
A 7
1.0 3.0 1.0
A 8
1.0 1.0 3.0
One species - starting as diamond cubic lattice
#volume 64.0
temperature 1.0
cutoff 1.0
timestep 0.01
steps 1000
traj 0 10 0
#traj 500 100 0
stats every 100
stack size 100
print every 100
job time 7200.0
close time 100.0
ensemble nvt mdvv
conf zero
nfold 2,2,2
One species - starting as diamond cubic lattice
A 1.0 0.0 8 0
A A dpd 25.0 1.0 4.0
# Local ordering for beads of species: A
# dimx, dimy, dimz= 8.0000000000000000 8.0000000000000000 8.0000000000000000
# snapshot number, q, sk
1 1.000000E+00 1.000000E+00
2 9.655988E-01 9.983900E-01
3 8.328311E-01 9.943825E-01
4 5.576891E-01 9.920264E-01
5 3.822875E-01 9.908347E-01
6 3.333387E-01 9.899954E-01
7 2.924687E-01 9.884140E-01
8 2.088436E-01 9.867583E-01
9 1.918481E-01 9.849730E-01
10 1.342873E-01 9.830251E-01
11 1.767485E-01 9.835370E-01
12 1.892161E-01 9.854308E-01
13 1.749140E-01 9.872711E-01
14 1.490145E-01 9.880391E-01
15 1.523878E-01 9.867670E-01
16 1.053237E-01 9.843228E-01
17 1.026490E-01 9.827817E-01
18 1.057525E-01 9.832584E-01
19 1.626147E-01 9.851831E-01
20 1.629386E-01 9.882356E-01
21 2.019539E-01 9.904825E-01
22 2.362456E-01 9.914645E-01
23 2.316278E-01 9.906415E-01
24 1.812034E-01 9.882595E-01
25 6.242488E-02 9.862133E-01
26 2.427986E-02 9.857918E-01
27 6.345471E-02 9.853107E-01
28 1.205525E-03 9.854216E-01
29 1.341872E-03 9.848937E-01
30 2.630008E-02 9.839716E-01
31 -3.033146E-02 9.833301E-01
32 6.626584E-03 9.826553E-01
33 2.761218E-02 9.830598E-01
34 4.759573E-02 9.845413E-01
35 5.307131E-02 9.864393E-01
36 5.445786E-02 9.867548E-01
37 1.711075E-01 9.870574E-01
38 1.265103E-01 9.870730E-01
39 9.441359E-02 9.866404E-01
40 7.403628E-02 9.871061E-01
41 8.127528E-02 9.887956E-01
42 1.109629E-01 9.898710E-01
43 1.609615E-01 9.899626E-01
44 1.737683E-01 9.896299E-01
45 1.335373E-01 9.887266E-01
46 1.519291E-01 9.875421E-01
47 1.496011E-01 9.864387E-01
48 1.030202E-01 9.861172E-01
49 6.936915E-02 9.859001E-01
50 7.178596E-03 9.850712E-01
51 -6.285664E-02 9.839900E-01
52 -9.091124E-03 9.839400E-01
53 4.894677E-02 9.853947E-01
54 7.888128E-02 9.876096E-01
55 8.285906E-02 9.884312E-01
56 8.505728E-02 9.868824E-01
57 6.781380E-02 9.847174E-01
58 5.306247E-02 9.839093E-01
59 5.050537E-02 9.849950E-01
60 7.860809E-02 9.865393E-01
61 9.849026E-02 9.881854E-01
62 1.066716E-01 9.891651E-01
63 8.805740E-02 9.884468E-01
64 1.045904E-01 9.871158E-01
65 8.977694E-02 9.852848E-01
66 5.350253E-02 9.840631E-01
67 2.858617E-02 9.842586E-01
68 4.736995E-02 9.856341E-01
69 8.104986E-02 9.864424E-01
70 1.020068E-01 9.871744E-01
71 1.170797E-01 9.879326E-01
72 1.365438E-01 9.886645E-01
73 1.382472E-01 9.890716E-01
74 1.149047E-01 9.888346E-01
75 6.268460E-02 9.882175E-01
76 4.625669E-02 9.869271E-01
77 1.633255E-02 9.858257E-01
78 7.265452E-02 9.852137E-01
79 5.460076E-02 9.852577E-01
80 6.783211E-02 9.857661E-01
81 4.952105E-02 9.852488E-01
82 6.424486E-02 9.851967E-01
83 6.248702E-02 9.866739E-01
84 8.541328E-02 9.883806E-01
85 1.341193E-01 9.893769E-01
86 1.830283E-01 9.895427E-01
87 2.099508E-01 9.894556E-01
88 1.684130E-01 9.889151E-01
89 1.506723E-01 9.880698E-01
90 9.588703E-02 9.873247E-01
91 6.160726E-02 9.863405E-01
92 5.526196E-02 9.856241E-01
93 6.741701E-02 9.854853E-01
94 -1.387784E-02 9.859535E-01
95 -2.738465E-03 9.858982E-01
96 -2.515220E-02 9.858201E-01
97 3.517364E-02 9.863703E-01
98 -5.866022E-03 9.862202E-01
99 3.245957E-02 9.854355E-01
100 1.102134E-01 9.853250E-01
101 1.279534E-01 9.861740E-01
# <q> = 0.127967E+00
# error = 0.165168E-01
# <s_k> = 0.986985E+00
# error = 0.277033E-03
Number of nodes used in calculations?
Which species number has to be analyzed?
total number of beads: 64
number of beads by species: 64
number of analyzed beads: 64
<q> = 0.127967E+00
error = 0.165168E-01
<s_k> = 0.986985E+00
error = 0.277033E-03
.. _dlmeso_tetrahedral:
WIP: Analysis of local tetrahedral ordering for DL_MESO_DPD
.. sidebar:: Software Technical Information
Documentation Tool
Application Documentation
See the Source Code section
Relevant Training Material
See the Testing section
.. contents:: :local:
Purpose of Module
This module, ``tetrahedral.f90``, is a postprocessing utility for DL_MESO_DPD,
the Dissipative Particle Dynamics (DPD) code from the DL_MESO_ package.
It processes the trajectory (HISTORY) files and analyzes the local tetrahedral
ordering, a feature that is relevant, for example, in water-like systems.
The local ordering in liquid water can be assessed considering the coordinates
of oxygen atoms [Duboue2015]_. In particular, for each oxygen, its four nearest neighbouring
oxygens are considered, whereas the hydrogens are disregarded.
At the mesoscale level, the user will select one (appropriate) bead species and analyze its local
Given a particle :math:`j`, we first find its *four nearest neighbours*
(n.n.). Then, an *orientational tetrahedral order parameter* is built using
where :math:`i,k` are n.n. of :math:`j` and :math:`\psi_{ik}=\theta_{ijk}` is the angle [1]_
between the particles :math:`i`, :math:`j` and :math:`k`.
Of course, the quantity is then averaged over the central particle :math:`j` and over time.
A *translational tetrahedral order parameter*, :math:`S_k`, is defined as
:math:`S_k=1-\frac{1}{3}\sum_{i=1}^4 \frac{(r_i - \bar{r})^2}{4\bar{r}^2}`,
where :math:`i` is a n.n. of :math:`j` and :math:`\bar{r}=\frac{1}{4}\sum_{i=1}^4 r_i`.
Concerning the limiting values of these parameters:
in a regular tetrahedron (if the four vertices are referred to the
center of the solid) one has :math:`q=S_k=1`.
In an ideal gas, where the angle :math:`\psi_{ik}` is randomly distributed,
:math:`q\simeq0`. On the other side, :math:`S_k\simeq 0` if the density
fluctuations are large enough.
As a result of the analysis, a file `TETRADAT` is produced, whose columns are
:math:`\textrm{snapshot index}, q, S_k`, the instantaneous values of the order
parameters defined above. At the end of the file, the averages and
standard errors (computed assuming the snapshots are uncorrelated) of both order parameters are given.
Background Information
The base code for this module is DL_MESO_DPD, the Dissipative Particle
Dynamics code from the mesoscopic simulation package DL_MESO_,
developed by M. Seaton at Daresbury Laboratory.
This open source code is available from STFC under both academic (free) and
commercial (paid) licenses. The module is to be used with DL_MESO_
in its second to last released version, version 2.6 (dating November 2015).
A variant of this module to be used with its last released version,
version 2.7 (released December 2018), will be provided soon.
The utility ``tetrahedral.f90`` is compiled with the available
Fortran90 compiler [2]_, e.g.:
``gfortran -o tetrahedral.exe tetrahedral.f90``
and the executable must be in the same directory of the HISTORY file.
The user is asked to provide the number of nodes used to run the
simulation and the number of the species for which ordering has to be
analyzed. To input the user-defined parameters one can enter them from the
keyboard or write them
into a text file (say, input.txt), one per line
(in the right order) and run the program in this way:
``tetrahedral.exe < input.txt``
Below we propose a test where a fluid is prepared in a ordered configuration
(diamond cubic lattice)
and rapidly goes into an orientationally disordered one.
Run the DL_MESO_DPD simulation on a single node (serial run)
using for the CONTROL file
.. literalinclude:: ./CONTROL
for the FIELD file
.. literalinclude:: ./FIELD
and for the CONFIG file
.. literalinclude:: ./CONFIG
This configuration corresponds to a diamond cubic lattice [3]_.
Analyzing the trajectory (HISTORY) file with ``tetrahedral.exe`` and input :math:`1`
for both requirements, this output is printed on the standard output
.. literalinclude:: ./out
The output file ``TETRADAT``
.. literalinclude:: ./TETRADAT
:lines: 1-13
containing the values of :math:`q` and :math:`S_k` for each snapshot and their
averages is produced too.
One can see that in the initial snapshot both order parameters detect an
ordered state (i.e., :math:`S_k=q=1`).
With the evolution in time, since the system is a dilute fluid,
the orientational ordering is rapidly lost (i.e., :math:`q\simeq 0`).
On the other side, the translational order parameter stays close to one since the density of the
system is roughly uniform.
Source Code
.. literalinclude:: ./tetrahedral.f90
:language: fortran
.. Here are the URL and references used
.. _DL_MESO:
.. [Duboue2015] E. Duboué-Dijon, A. Laage, *Characterization of the
local structure in liquid water by various order parameters*,
J. Phys. Chem. B, **119**, 8406 (2015).
.. [1] The angle
where :math:`\vec{r_{ij}} = \vec{r_i} -\vec{r_j}` and :math:`r=|\vec{r}|`.
.. [2] Compilation has been tested with the GNU compiler GCC, version 8.2.1.
.. [3] The diamond cubic crystal lattice ( is a repeating pattern of
8 atoms. Their coordinates may be given as:
:math:`A=(0,0,0)`, :math:`B=(0,2,2)`, :math:`C=(2,0,2)`, :math:`D=(2,2,0)`, :math:`E=(3,3,3)`,
:math:`F=(3,1,1)`, :math:`G=(1,3,1)`, and :math:`H=(1,1,3)` in a unit cubic cell of side :math:`L=4`.
One can check that, with the minimum image convention, each
particle has its 4 closest neighbours at a distance :math:`\sqrt{3}`, and all the angles are
For this configuration (also if repeated periodically along the three Cartesian axis), :math:`q=1` and :math:`S_k=1`.
PROGRAM tetrahedral
! module to analyze tetrahedral ordering in dl_meso HISTORY files
! authors - m. a. seaton & s. chiacchiera, january 2018 (tidied up on march 2019)
REAL(KIND=dp), PARAMETER :: pi=3.141592653589793_dp
CHARACTER(80) :: text, a2
CHARACTER(8), ALLOCATABLE :: namspe (:), nammol (:)
CHARACTER(6) :: chan
CHARACTER(8) :: a1
INTEGER, ALLOCATABLE :: ltp (:), ltm (:), beads (:), bonds (:), nspec (:)
INTEGER :: nrtout
INTEGER :: chain, imol, ibead, ioerror, i, numtraj, j, k, l, m, nmoldef, ibond
INTEGER :: nspe, numnodes, nbeads, nusyst, nmbeads, nsyst, nbonds, numbond, global, species, molecule
INTEGER :: nummol, lfrzn, rnmol, keytrj, srfx, srfy, srfz
INTEGER :: indx, nav
INTEGER :: n1, n2, n3, n4
INTEGER :: bead1, bead2
REAL(KIND=dp), ALLOCATABLE :: xxx (:), yyy (:), zzz (:)
REAL(KIND=dp), ALLOCATABLE :: nmol (:)
REAL(KIND=dp) :: volm, dimx, dimy, dimz, shrdx, shrdy, shrdz
REAL(KIND=dp) :: amass, rcii
REAL(KIND=dp) :: time, mbeads, mglobal, x, y, z, vx, vy, vz, fx, fy, fz, rsq
REAL(KIND=dp) :: r1, r2, r3, r4
LOGICAL :: eof
! Variables for tetrahedral ordering
INTEGER :: nnlab(4), npart, sp, count
REAL(KIND=dp) :: qtetra, stetra
REAL(KIND=dp) :: q, sk, q_sum, sk_sum, q_ave, sk_ave
REAL(KIND=dp) :: q2_sum, sk2_sum, q2_ave, sk2_ave
! Get number of nodes
WRITE (*,*) "Number of nodes used in calculations?"
READ (*,*) numnodes
ALLOCATE (beads (numnodes), bonds (numnodes))
! Determine if HISTORY files exist
IF (numnodes>1) THEN
INQUIRE (file = 'HISTORY000000', EXIST = eof)
INQUIRE (file = 'HISTORY', EXIST = eof)
IF (.NOT. eof) THEN
WRITE (*,*) "ERROR: cannot find HISTORY files"
! First reading, where the number of beads, molecules and bonds are determined
! Arrays are filled with names of particles and molecules
! If multiple HISTORY files are present, it is checked they are compatible
numbond = 0
DO j = 1, numnodes
WRITE (chan, '(i6.6)') j-1
IF (numnodes>1) THEN
OPEN (ntraj+j-1, file = 'HISTORY'//chan, access = 'sequential', form = 'unformatted', status = 'unknown')
OPEN (ntraj, file = 'HISTORY', access = 'sequential', form = 'unformatted', status = 'unknown')
IF (j == 1) THEN
READ (ntraj+j-1) nspe, nmoldef, nusyst, nsyst, nbeads, nbonds
READ (ntraj+j-1) dimx, dimy, dimz, volm
READ (ntraj+j-1) keytrj, srfx, srfy, srfz
READ (ntraj+j-1) n1, n2, n3, n4, nbeads, nbonds
READ (ntraj+j-1) r1, r2, r3, r4
IF (n1 /= nspe .OR. n2 /= nmoldef .OR. n3 /= nusyst .OR. n4 /= nsyst &
.OR. r1 /= dimx .OR. r2 /= dimy .OR. r3 /= dimz .OR. r4 /= volm) THEN
WRITE (*,*) "ERROR: HISTORY files do not refer to the same system!"
READ (ntraj+j-1) n1, n2, n3, n4
IF (n1 /= keytrj .OR. n2 /= srfx .OR. n3 /= srfy .OR. n4 /= srfz) THEN
WRITE (*,*) "ERROR: HISTORY files do not refer to the same system!"
beads (j) = nbeads
bonds (j) = nbonds
numbond = numbond + nbonds
END DO ! loop over nodes
! IF (srfx == 3 .OR. srfy == 3 .OR. srfz == 3) THEN
! WRITE (*,*) "ERROR: System under shear, not implemented yet!"
ALLOCATE (namspe (nspe), nammol (nmoldef), nspec (nspe)) !NB: nspec here counts ALL beads of a type, not only unbonded ones
ALLOCATE (xxx (1:nsyst), yyy (1:nsyst), zzz (1:nsyst))
ALLOCATE (ltp (1:nsyst))
DO j = 1, numnodes
DO i = 1, nspe
IF (j == 1) THEN
READ (ntraj+j-1) namspe (i), amass, rcii, lfrzn
READ (ntraj+j-1) a1, amass, rcii, lfrzn
IF (a1 /= namspe (i))THEN
WRITE (*,*) "ERROR: HISTORY files do not refer to the same system!"
IF (nmoldef>0) THEN
DO i = 1, nmoldef
IF (j==1) THEN
READ (ntraj+j-1) nammol (i)
READ (ntraj+j-1) a1
IF (a1 /= nammol (i))THEN
WRITE (*,*) "ERROR: HISTORY files do not refer to the same system!"
IF (j == 1) THEN
READ (ntraj+j-1) text
READ (ntraj+j-1) a2
IF (a2 /= text) THEN
WRITE (*,*) "ERROR: HISTORY files do not refer to the same system!"
ENDDO ! end of loop over nodes
DO j = 1, numnodes
CLOSE (ntraj+j-1)
! Second reading, where arrays are filled with properties of beads and molecules.
! Then, the snapshots of trajectories are read.
DO j = 1, numnodes
WRITE (chan, '(i6.6)') j-1
IF (numnodes>1) THEN
OPEN (ntraj+j-1, file = 'HISTORY'//chan, access = 'sequential', form = 'unformatted', status = 'unknown')
OPEN (ntraj, file = 'HISTORY', access = 'sequential', form = 'unformatted', status = 'unknown')
READ (ntraj+j-1) !nspe, nmoldef, nusyst, nsyst, nbeads, nbonds
READ (ntraj+j-1) !dimx, dimy, dimz, volm
READ (ntraj+j-1) !keytrj, srfx, srfy, srfz
DO i = 1, nspe
READ (ntraj+j-1) !namspe (i), amass, rcii, lfrzn
DO i = 1, nmoldef
READ (ntraj+j-1) !nammol (i)
READ (ntraj+j-1) !text
ibond = 0 !counter for bonds
nspec (:) = 0 ! populations
! fill in arrays for beads and bonds
DO j = 1, numnodes
!Build ltp, ltm, mole
DO i = 1, beads (j)
READ (ntraj+j-1) global, species, molecule, chain
ltp (global) = species
nspec (species) = nspec (species) + 1
IF (bonds (j)>0) THEN
DO i = 1, bonds (j)
ibond = ibond + 1
READ (ntraj+j-1) bead1, bead2
END DO ! over nodes
IF (ibond /= numbond) THEN
WRITE (*,*) "ERROR: mismatch in number of bonds!"
! Find number of beads for which trajectories are needed
WRITE(*,*) "Which species number has to be analyzed?"
READ(*,*) sp
IF (sp<0 .OR. sp>nspe) THEN
WRITE(*,*) "error: undefined species!"
npart = nspec (sp)
WRITE(*,*) "total number of beads: ", nsyst
WRITE(*,*) "number of beads by species: ", nspec
WRITE(*,*) "number of analyzed beads: ", npart