Commit 5efba403 authored by Jony's avatar Jony

Added DL_MESO_DPD folder

parent 1bb87982
Two kind of molecules, branched and dimer
0 1
3.0 0.0 0.0
0.0 3.0 0.0
0.0 0.0 3.0
B 1
0.0 0.0 0.0
A 2
0.0 0.2 0.0
C 3
0.0 0.4 0.0
A 4
0.2 0.2 0.0
B 5
0.0 0.0 0.3
D 6
0.0 0.0 0.1
Two kinds of molecules: branched and dimer
volume 3.0 3.0 3.0
temperature 1.0
cutoff 1.0
timestep 0.01
steps 1000
equilibration steps 0
traj 0 100 0
stats every 100
stack size 100
print every 100
job time 1000.0
close time 10.0
ensemble nvt mdvv
conf origin zero
ewald sum 1.0 5 5 5
bjerrum 1.0
smear gauss equal
finish
Two kinds of molecules: branched and dimer
SPECIES 4
A 1.0 0.2 0 0
B 1.0 -1.0 0 0
C 1.0 0.6 0 0
D 1.0 1.0 0 0
MOLECULES 2
BRANCH
nummols 10
beads 4
B 0.0 0.0 0.0
A 0.0 0.2 0.0
C 0.0 0.4 0.0
A 0.2 0.2 0.0
bonds 3
harm 1 2 5.0 0.25
harm 2 3 5.0 0.25
harm 2 4 5.0 0.25
finish
BD
nummols 10
beads 2
B 0.0 0.0 0.3
D 0.0 0.0 0.1
bonds 1
harm 1 2 5.0 0.25
finish
INTERACTIONS 4
A A dpd 25.0 1.0 4.5
B B dpd 25.0 1.0 4.5
C C dpd 25.0 1.0 4.5
D D dpd 25.0 1.0 4.5
CLOSE
1 5.610985E-02 1.656239E-01 -1.397198E-01 4.000000E-02 1.855601E-03
2 4.642350E-01 -9.118123E-01 9.524455E-02 1.599716E+00 3.911064E-02
3 -2.027409E+00 -1.131613E+00 3.758015E-01 4.747804E-01 2.048949E-01
4 1.319446E+00 1.476795E-01 1.063918E+00 8.070183E-01 1.072100E-01
5 2.010414E+00 -4.817778E-01 2.350176E+00 9.949271E-01 3.628594E-01
6 -1.451075E+00 9.994575E-01 3.334998E+00 1.227528E+00 5.269164E-01
7 -2.657225E-01 -7.582258E-01 2.300250E-01 7.640699E-01 2.586764E-02
8 5.496076E-01 -8.682396E-01 1.023472E+00 8.465722E-01 7.790383E-02
9 1.088395E+00 -2.984139E-01 -2.146073E+00 9.144278E-01 2.177513E-01
10 -8.548484E-01 -7.038363E-01 9.717777E-01 7.643206E-01 8.038901E-02
11 2.299682E+00 -3.035233E+00 -2.349568E+00 8.803798E-01 7.415426E-01
1 0.000000E+00 0.000000E+00 -2.000000E-01 4.000000E-02 1.481481E-03
1 2.869309E-02 -4.152495E-01 1.999305E-01 1.040000E-01 7.897321E-03
2 -1.875110E+00 -9.240273E-01 -1.062192E+00 1.335149E+00 2.036340E-01
3 1.333150E+00 2.371575E+00 4.829357E-02 1.537261E+00 2.742218E-01
4 -7.447839E-02 7.132180E-01 -5.185719E-01 1.806593E+00 2.900532E-02
5 8.163357E+00 -2.563583E+00 -2.268061E-01 1.699112E+00 2.713474E+00
6 1.909976E+00 -2.072637E+00 1.328292E+00 1.408593E+00 3.595626E-01
7 -1.606210E+00 -3.637945E+00 1.086719E+00 1.279568E+00 6.294633E-01
8 -7.687367E-01 -6.351486E+00 1.252305E+00 1.355228E+00 1.574096E+00
9 1.428639E-01 -3.620968E+00 3.834862E-01 1.132792E+00 4.918104E-01
10 2.524788E+00 1.276311E+00 1.641325E+00 1.187329E+00 3.962026E-01
11 1.267564E+00 -5.540963E-01 -2.349736E+00 1.672068E+00 2.753703E-01
1 4.000000E-02 3.200000E-01 0.000000E+00 1.040000E-01 3.851852E-03
This diff is collapsed.
.. _example:
################################################
Analysis of charge dipole moments in DL_MESO_DPD
################################################
.. sidebar:: Software Technical Information
Language
FORTRAN 90
Licence
Documentation Tool
Latex-generated .pdf file
Application Documentation
Click to see the manual_ with more details
Relevant Training Material
.. Add a link to any relevant training material.
.. contents:: :local:
Purpose of Module
_________________
This module, ``gen_dipole.f90``, is a generalization of the ``dipole.f90`` post-processing
utility of DL_MESO_DPD, the Dissipative Particle Dynamics (DPD) code from the DL_MESO_ package.
It processes the trajectory (HISTORY) files to obtain the charge dipole moments
of all the (neutral) molecules in the system.
It produces files `dipole_*` containing the time evolution of relevant
quantities (see below). In the case of a single molecular species, it also prints
to the standard output the Kirkwood number :math:`g_k` and the relative electric
permittivity :math:`\epsilon_r` for this species, together with an estimate for their errors (standard error).
The module can be applied to systems including molecules with a generic charge structure, as long
as each molecule is neutral (otherwise the charge dipole moment would be frame-dependent).
The charge dipole moment of a neutral molecule is :math:`\vec{p}_{mol}=\sum_{i\in mol}q_i \vec{r}_i` where
:math:`\vec{r}_i` are the bead positions and :math:`q_i` their charges. The
total charge dipole moment of the simulated volume :math:`V` is
:math:`\vec{P}=\sum_{mol\in V} \vec{p}_{mol}`.
If more than one molecular species are present, one can split :math:`\vec{P}` into the different species contributions.
In general:
For any molecular species a file `dipole_{molecule name}` is produced, whose columns are
:math:`\textrm{snapshot index}, P_x, P_y, P_z, \sum_{i=1}^{N_{mol}}\frac{\vec{p}_i ^2}{N_{mol}},\frac{\vec{P} ^2}{V}`.
It is intended that for any quantity the contribution given *from the
species {molecule name}* is reported (i.e., the sums are restricted to
molecules of a single type).
Possible uses of the output files are: monitoring the polarization in response to an
external electric field, measuring the fluctuations in molecular/total charge
dipole moments.
Extra output for a single molecular species:
The Kirkwood number for a pure system is
:math:`g_k=\frac{\langle\vec{P}^2\rangle}{N_{mol}\langle \vec{p}^2\rangle}`,
where :math:`\langle\dots\rangle` indicates an average over trajectories.
If the dipoles' orientations are not correlated, then :math:`g_k\simeq 1`.
Also, an estimate for the relative dielectric permittivity of the medium is obtained from linear response
theory: :math:`\epsilon_r= 1 + \frac{4 \pi}{3} l_B \frac{\langle\vec{P}^2\rangle}{V}`,
where :math:`l_B` is Bjerrum length.
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 last released version, version 2.6 (dating November 2015).
Testing
_______
The present module ``gen_dipole.f90`` is compiled with the available Fortran90 compiler, e.g.:
``gfortran -o gen_dipole.exe gen_dipole.f90``
and the executable must be in the same directory of the HISTORY* files to be
analyzed. The user is asked to provide the number of nodes used to run the
simulation, the electric charges for all the types of beads
and the Bjerrum length.
As a test, we suggest two examples involving two (toy) molecular species: a
branched one (four beads, T-shaped) and a simple dimer. All the beads carry charges.
In the first case 10 molecules of each type are present and are followed for a few time steps.
In the second case it is suggested to analyze a single snapshot with just two
molecules and all the beads sitting at user-defined positions (via the CONFIG
file).
Four type of beads are used with charges :math:`q_A=0.2, q_B=-1, q_C=0.6,q_D=1`;
the Bjerrum length is fixed as :math:`l_B=1`.
The bonding connections in the two molecules are pictorially given below:
::
B - A - C B - D
|
A
**First test**
Run the DL_MESO_DPD simulation on a single node (serial run)
using for the CONTROL file
.. literalinclude:: ./CONTROL
and for the FIELD file
.. literalinclude:: ./FIELD
Analyzing the HISTORY file with `gen_dipole.exe`, this output is printed on the standard output
.. literalinclude:: ./out-d
The first line shows the histogram of cluster sizes: in this case,
it correctly gives 10 molecules of two beads, and 10 molecules of 4 beads.
Since internally the module checks that each molecule is a connected cluster [1]_,
this line should always give a histogram with the molecule sizes
(by default, shown up to ten beads).
The `dipole_BD` file is
.. literalinclude:: ./dipole_BD
and the `dipole_BRANCH` one is
.. literalinclude:: ./dipole_BRANCH
If instead the simulation is run on multiple nodes, only the results for the
first snapshot will be unchanged (i.e., first line of each `dipole_*` file),
the other results will vary because a different sequence of random numbers
will enter the time evolution of the system.
**Second test**
Run DL_MESO_DPD using the same CONTROL and FIELD files as above, with the only changes:
- `"steps 1000"` changes into `"steps 1"`
- `"nummols 10"` changes into `"nummols 1"` (NB: appears twice)
Also, use this CONFIG file that will initially put the molecules
branches aligned with the cartesian axes
.. literalinclude:: ./CONFIG
where the identity of each bead is fixed by the FIELD file and is shown below
::
B(1) - A(2) - C(3) B(5) - D(6)
|
A(4)
One can easily check that the dipole of each molecule is as expected:
:math:`\vec{p}_{BRANCH}=(0.04,0.32,0), \quad \vec{p}_{BD}=(0,0,-0.2)~.`
In fact: the `dipole_BD` file is
.. literalinclude:: ./dipole_BD-2
and the `dipole_BRANCH` one is
.. literalinclude:: ./dipole_BRANCH-2
The results of this test do not depend on the number of nodes used to run the simulation.
Source Code
___________
.. literalinclude:: ./gen_dipole.f90
:language: fortran
:linenos:
.. _manual: ./source/mandip.pdf
.. Here are the URL references used
.. _DL_MESO: http://www.scd.stfc.ac.uk/SCD/support/40694.aspx
.. _ReST: http://docutils.sourceforge.net/docs/user/rst/quickref.html
.. [1] Disambiguation on the concept of molecule. In DL\_MESO a *defined molecule*
is a set of beads, which can be bonded or not.
For the purpose of this module it is *required* that each molecule is a
connected cluster (via stretching bonds).
In fact, this, together with the reasonable assumption that each stretching
bond cannot be stretched to more than half the system linear size, allows
to univocally define the charge dipole moment of each molecule.
.. www.ccp5.ac.uk/DL_MESO
Simple test
volume 3.0 3.0 3.0
temperature 1.0
cutoff 1.0
timestep 0.01
steps 6
equilibration steps 2
traj 2 2 0
stats every 2
stack size 2
print every 2
job time 100.0
close time 10.0
#surface shear y
#surface frozen x
#surface hard x
ensemble nvt mdvv
finish
Simple test
SPECIES 3
A 1.0 0.0 1 0
B 1.0 0.0 0 0
C 1.0 0.0 0 0
MOLECULES 2
AB
nummols 1
beads 2
A 0.0 0.0 0.0
B 0.1 0.0 0.0
bonds 1
harm 1 2 5.0 0.0
finish
AC
nummols 1
beads 2
A 0.0 0.0 0.0
C 0.1 0.0 0.0
bonds 1
harm 1 2 3.0 0.0
finish
INTERACTIONS 3
A A dpd 25.0 1.0 4.5
B B dpd 25.0 1.0 4.5
C C dpd 25.0 1.0 4.5
#EXTERNAL
#shear 3.0 0.0 0.0
CLOSE
PROGRAM format_history
!***********************************************************************************
!
! module to format dl_meso HISTORY files
!
! authors - m. a. seaton & s. chiacchiera, february 2017
!
!**********************************************************************************
IMPLICIT none
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND (15, 307)
INTEGER, PARAMETER :: ntraj=10,nuser=5
CHARACTER(80) :: text, a2
CHARACTER(8), ALLOCATABLE :: namspe (:), nammol (:)
CHARACTER(6) :: chan
CHARACTER(8) :: a1
INTEGER, ALLOCATABLE :: ltp (:), ltm (:), mole (:), beads (:), bonds (:), bndtbl (:,:)
INTEGER, ALLOCATABLE :: nbdmol (:), nbomol (:)
INTEGER :: chain, imol, ioerror, i, k, j, nmoldef, ibond
INTEGER :: nspe, nbeads, nusyst, nsyst, nbonds, global, species, molecule, numnodes, numbond
INTEGER :: nummol, lfrzn, rnmol, keytrj, srfx, srfy, srfz
INTEGER :: bead1, bead2
INTEGER :: n1, n2, n3, n4
INTEGER :: nform
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
REAL(KIND=dp) :: r1, r2, r3, r4
LOGICAL :: eof, lcomm, lmcheck
! Switches for commenting and checking molecules
lcomm = .TRUE.
lmcheck = .TRUE.
! 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)
ELSE
INQUIRE (file = 'HISTORY', EXIST = eof)
END IF
IF (.NOT. eof) THEN
WRITE (*,*) "ERROR: cannot find HISTORY files"
STOP
END IF
! Open the output files
nform = ntraj + numnodes
DO j = 1, numnodes
WRITE (chan, '(i6.6)') j-1
IF (numnodes>1)THEN
OPEN (nform+j-1, file = 'HISTORY'//chan//"-F", status = 'replace')
ELSE
OPEN (nform+j-1, file = 'HISTORY'//"-F", status = 'replace')
END IF
END DO
! 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')
ELSE
OPEN (ntraj, file = 'HISTORY', access = 'sequential', form = 'unformatted', status = 'unknown')
END IF
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
ELSE
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!"
STOP
ENDIF
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!"
STOP
ENDIF
ENDIF
beads (j) = nbeads
bonds (j) = nbonds
numbond = numbond + nbonds
IF (lcomm) WRITE (nform+j-1,*) "# nspe, nmoldef, nusyst, nsyst, nbeads, nbonds"
WRITE (nform+j-1,*) nspe, nmoldef, nusyst, nsyst, nbeads, nbonds
IF (lcomm) WRITE (nform+j-1,*) "# dimx, dimy, dimz, volm"
WRITE (nform+j-1,97) dimx, dimy, dimz, volm
IF (lcomm) WRITE (nform+j-1,*) "# keytrj, srfx, srfy, srfz"
WRITE (nform+j-1,*) keytrj, srfx, srfy, srfz
END DO ! loop over nodes
ALLOCATE (namspe (nspe), nammol (nmoldef))
IF (lmcheck) THEN
ALLOCATE (ltp (1:nsyst), ltm (1:nsyst), mole (1:nsyst))
ALLOCATE (nmol (1:nmoldef), nbdmol (1:nmoldef), nbomol (1:nmoldef))
ALLOCATE (bndtbl (numbond, 2))
ENDIF
DO j = 1, numnodes
IF (lcomm) WRITE (nform+j-1,*) "# SPECIES:"
IF (lcomm) WRITE (nform+j-1,*) "# namspe, amass, rcii, lfrzn"
DO i = 1, nspe
IF (j == 1) THEN
READ (ntraj+j-1) namspe (i), amass, rcii, lfrzn
ELSE
READ (ntraj+j-1) a1, amass, rcii, lfrzn
IF (a1 /= namspe (i))THEN
WRITE (*,*) "ERROR: HISTORY files do not refer to the same system!"
STOP
ENDIF
ENDIF
WRITE (nform+j-1,96) namspe (i), amass, rcii, lfrzn
END DO
IF (nmoldef>0) THEN
IF (lcomm) WRITE (nform+j-1,*) "# MOLECULES:"
IF (lcomm) WRITE (nform+j-1,*) "# nammol"
DO i = 1, nmoldef
IF (j==1) THEN
READ (ntraj+j-1) nammol (i)
ELSE
READ (ntraj+j-1) a1
IF (a1 /= nammol (i))THEN
WRITE (*,*) "ERROR: HISTORY files do not refer to the same system!"
STOP
ENDIF
END IF
WRITE (nform+j-1,*) nammol (i)
END DO
END IF
IF (j == 1) THEN
READ (ntraj+j-1) text
ELSE
READ (ntraj+j-1) a2
IF (a2 /= text) THEN
WRITE (*,*) "ERROR: HISTORY files do not refer to the same system!"
STOP
ENDIF
ENDIF
IF (lcomm) WRITE (nform+j-1,*) "# Simulation name:"
WRITE (nform+j-1,*) text
ENDDO ! end of loop over nodes
DO j = 1, numnodes
CLOSE (ntraj+j-1)
END DO
! Second reading, where (if required) 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')
ELSE
OPEN (ntraj, file = 'HISTORY', access = 'sequential', form = 'unformatted', status = 'unknown')
END IF
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
END DO
DO i = 1, nmoldef
READ (ntraj+j-1) !nammol (i)
END DO
READ (ntraj+j-1) !text
END DO
nummol = 0 !counter for number of molecules
ibond = 0 !counter for bonds
! fill in arrays for beads and bonds
DO j = 1, numnodes
IF (lcomm) WRITE (nform+j-1,*) "# BEADS:"
IF (lcomm) WRITE (nform+j-1,*) "# global, species, molecule, chain"
IF (lmcheck) THEN
!Build ltp, ltm, mole
DO i = 1, beads (j)
READ (ntraj+j-1) global, species, molecule, chain
ltp (global) = species
ltm (global) = molecule
mole (global) = chain
nummol = MAX (nummol, chain)
WRITE (nform+j-1,*) global, species, molecule, chain
END DO
ELSE
DO i = 1, beads (j)
READ (ntraj+j-1) global, species, molecule, chain
WRITE (nform+j-1,*) global, species, molecule, chain
END DO
ENDIF
IF (bonds (j)>0) THEN
IF (lcomm) WRITE (nform+j-1,*) "# BONDS:"
IF (lcomm) WRITE (nform+j-1,*) "# extremes of the bond"
IF (lmcheck) THEN
! Build bndtbl
DO i = 1, bonds (j)
ibond = ibond + 1
READ (ntraj+j-1) bead1, bead2
bndtbl (ibond, 1) = bead1
bndtbl (ibond, 2) = bead2
WRITE (nform+j-1,*) bead1, bead2
END DO
ELSE
DO i = 1, bonds (j)
READ (ntraj+j-1) bead1, bead2
WRITE (nform+j-1,*) bead1, bead2
END DO
END IF
END IF
</