Commit 2b4a5b4c authored by Jony's avatar Jony

Removed previous folders

parent 5efba403
.. _example:
####################
[Wip] 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
.. Provide a link to any documentation.
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 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).
.. www.ccp5.ac.uk/DL_MESO
.. M. A. Seaton, R. L. Anderson, S. Metz and W. Smith, Mol. Sim. 39 (10), 796 (2013).
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. To test ...
Source Code
___________
.. Here link the source code that was created for the module.
.. 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
.. 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 - s. chiacchiera & m. a. seaton, february 2017
!
!**********************************************************************************
IMPLICIT none
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND (15, 307)
INTEGER, PARAMETER :: ntraj=10,nuser=5!,nrtout=9
CHARACTER(80) :: text
CHARACTER(8), ALLOCATABLE :: namspe (:), nammol (:)
! CHARACTER(6) :: chan
INTEGER, ALLOCATABLE :: ltp (:), ltm (:), mole (:)!, bndtbl (:,:), beads (:), bonds (:)
INTEGER, ALLOCATABLE :: nbdmol (:)!, nbond (:)
INTEGER :: chain, imol, ioerror, i, k, nmoldef!, ntop, numtraj, j, l, m, ibond, ibead, mxgrid
INTEGER :: nspe, nbeads, nusyst, nsyst, nbonds, global, species, molecule !numnodes, numbonds
INTEGER :: nummol, lfrzn, rnmol, keytrj, srfx, srfy, srfz
INTEGER :: bead1, bead2
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
LOGICAL :: eof, lcomm, lmcheck
! Switches for commenting and checking molecules
lcomm = .TRUE.!.FALSE.!
lmcheck = .TRUE.!.FALSE.
! First reading, where the number of beads and molecules are determined
OPEN (ntraj, file = 'HISTORY', access = 'sequential', form = 'unformatted', status = 'unknown')
READ (ntraj) nspe, nmoldef, nusyst, nsyst, nbeads, nbonds
READ (ntraj) dimx, dimy, dimz, volm
READ (ntraj) keytrj, srfx, srfy, srfz
IF (lcomm) WRITE (*,*) "# nspe, nmoldef, nusyst, nsyst, nbeads, nbonds"
WRITE (*,*) nspe, nmoldef, nusyst, nsyst, nbeads, nbonds
IF (lcomm) WRITE (*,*) "# dimx, dimy, dimz, volm"
WRITE (*,97) dimx, dimy, dimz, volm
IF (lcomm) WRITE (*,*) "# keytrj, srfx, srfy, srfz"
WRITE (*,*) keytrj, srfx, srfy, srfz
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))
ENDIF
IF (lcomm) WRITE (*,*) "# SPECIES:"
IF (lcomm) WRITE (*,*) "# namspe, amass, rcii, lfrzn"
DO i = 1, nspe
READ (ntraj) namspe (i), amass, rcii, lfrzn
WRITE (*,96) namspe (i), amass, rcii, lfrzn
END DO
IF (nmoldef>0) THEN
IF (lcomm) WRITE (*,*) "# MOLECULES:"
IF (lcomm) WRITE (*,*) "# nammol"
DO i = 1, nmoldef
READ (ntraj) nammol (i)
WRITE (*,*) nammol (i)
END DO
END IF
READ (ntraj) text
IF (lcomm) WRITE (*,*) "# Simulation name:"
WRITE (*,*) text
IF (lcomm) WRITE (*,*) "# BEADS:"
IF (lcomm) WRITE (*,*) "# global, species, molecule, chain"
DO i = 1, nbeads
READ (ntraj) global, species, molecule, chain
WRITE (*,*) global, species, molecule, chain
END DO
IF (nbonds>0) THEN
IF (lcomm) WRITE (*,*) "# BONDS:"
IF (lcomm) WRITE (*,*) "# extremes of the bond"
DO i = 1, nbonds
READ (ntraj) bead1, bead2
WRITE (*,*) bead1, bead2
END DO
END IF
CLOSE (ntraj)
! Second reading, where (if required) arrays are filled with names and properties
! of beads and molecules. Then, the snapshots of trajectories are read.
OPEN (ntraj+nmoldef, file = 'HISTORY', access = 'sequential', form = 'unformatted', status = 'unknown')
READ (ntraj+nmoldef) !nspe, nmoldef, nusyst, nsyst, nbeads, nbonds
READ (ntraj+nmoldef) !dimx, dimy, dimz, volm
READ (ntraj+nmoldef) !keytrj, srfx, srfy, srfz
DO i = 1, nspe
READ (ntraj+nmoldef) !namspe (i), amass, rcii, lfrzn
END DO
DO i = 1, nmoldef
READ (ntraj+nmoldef) !nammol (i)
END DO
READ (ntraj+nmoldef) !text
! fill in arrays for beads
! ibond = 0
nummol = 0 !number of molecules
IF (lmcheck ) THEN
WRITE (*,*) "# Build and check ltp, ltm, mole"
DO i = 1, nbeads
READ (ntraj+nmoldef) global, species, molecule, chain
ltp (global) = species
ltm (global) = molecule
mole (global) = chain
nummol = MAX (nummol, chain)
WRITE(*,*) global, ltp (global), ltm (global), mole (global)
END DO
ELSE
DO i = 1, nbeads
READ (ntraj+nmoldef) !global, species, molecule, chain
END DO
ENDIF
IF (nbonds>0) THEN
DO i = 1, nbonds
READ (ntraj+nmoldef) !molecule, chain
END DO
END IF
IF (lmcheck) THEN
! !determine numbers of molecules and beads per molecule type ! I have removed the "bonds"
nmol = 0.0_dp
nbdmol = 0
! nbond = 0
chain = 0
DO i = 1, nbeads
IF (mole (i) /= chain) THEN
chain = mole (i)
imol = ltm (i)
nmol (imol) = nmol (imol) + 1.0_dp
END IF
IF (imol > 0) nbdmol (imol) = nbdmol (imol) + 1
END DO
DO i = 1, nmoldef
rnmol = NINT (nmol (i))
IF (rnmol>0) THEN
nbdmol (i) = nbdmol (i) / rnmol
! nbond (i) = nbond (i) / rnmol
END IF
END DO
!-
!Check of molecule beads and numbers
IF (nmoldef>0) THEN
WRITE (*,*) "# Check of molecules: nammol, ndbmol, nmol"
DO i = 1, nmoldef
WRITE (*,*) nammol (i), nbdmol (i), NINT(nmol(i))
END DO
WRITE (*,*) "# Total number of molecules = ",nummol
END IF
END IF
!reading trajectories
eof = .false.
k = 0
IF (lcomm) WRITE (*,*) "# --- TRAJECTORIES --- (key =", keytrj,")"
SELECT CASE (keytrj)
CASE (0)
IF (lcomm) WRITE (*,*) "# mglobal, x, y, z"
CASE(1)
IF (lcomm) WRITE (*,*) "# mglobal, x, y, z, vx, vy, vz"
CASE(2)
IF (lcomm) WRITE (*,*) "# mglobal, x, y, z, vx, vy, vz, fx, fy, fz"
END SELECT
DO WHILE (.true.)
READ (ntraj+nmoldef, IOSTAT=ioerror) time, mbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
IF (lcomm) WRITE (*,*) "# time, mbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz"
WRITE (*,98) time, mbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
IF (ioerror/=0) THEN
eof = .true.
IF (k==0) THEN
PRINT *, 'ERROR: cannot find trajectory data in HISTORY files'
STOP
END IF
EXIT
END IF
k = k + 1
IF (lcomm) WRITE (*,*) "# snapshot number:", k
nbeads = NINT (mbeads)
SELECT CASE (keytrj)
CASE (0)
DO i = 1, nbeads
READ (ntraj+nmoldef) mglobal, x, y, z
WRITE (*,99) mglobal, x, y, z
END DO
CASE (1)
DO i = 1, nbeads
READ (ntraj+nmoldef) mglobal, x, y, z, vx, vy, vz
WRITE (*,99) mglobal, x, y, z, vx, vy, vz
END DO
CASE (2)
DO i = 1, nbeads
READ (ntraj+nmoldef) mglobal, x, y, z, vx, vy, vz, fx, fy, fz
WRITE (*,99) mglobal, x, y, z, vx, vy, vz, fx, fy, fz
END DO
END SELECT
END DO
CLOSE (ntraj+nmoldef)
DEALLOCATE (namspe, nammol)
IF (lmcheck) DEALLOCATE (ltp, ltm, mole, nmol, nbdmol)
99 FORMAT(f10.1,2x,1p,9(e13.6,3x))
98 FORMAT(8(f10.3,3x))
97 FORMAT(4(f10.3,3x))
96 FORMAT(A9,3x,2(f10.3,3x),I2)
END PROGRAM format_history
# nspe, nmoldef, nusyst, nsyst, nbeads, nbonds
3 2 1 5 5 2
# dimx, dimy, dimz, volm
3.000 3.000 3.000 27.000
# keytrj, srfx, srfy, srfz
0 0 0 0
# SPECIES:
# namspe, amass, rcii, lfrzn
A 1.000 1.000 0
B 1.000 1.000 0
C 1.000 1.000 0
# MOLECULES:
# nammol
AB
AC
# Simulation name:
Simple test
# BEADS:
# global, species, molecule, chain
1 1 0 0
2 1 1 1
3 2 1 1
4 1 2 2
5 3 2 2
# BONDS:
# extremes of the bond
2 3
4 5
# Build and check ltp, ltm, mole
1 1 0 0
2 1 1 1
3 2 1 1
4 1 2 2
5 3 2 2
# Check of molecules: nammol, ndbmol, nmol
AB 2 1
AC 2 1
# Total number of molecules = 2
# --- TRAJECTORIES --- (key = 0 )
# mglobal, x, y, z
# time, mbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
0.000 5.000 3.000 3.000 3.000 0.000 0.000 0.000
# snapshot number: 1
1.0 1.471873E+00 1.525203E+00 1.507395E+00
2.0 1.364570E+00 2.228593E+00 2.475293E+00
3.0 1.300679E+00 2.256132E+00 2.340013E+00
4.0 2.306000E+00 2.535871E+00 2.718987E-01
5.0 2.292338E+00 2.576789E+00 2.972781E-01
# time, mbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
0.020 5.000 3.000 3.000 3.000 0.000 0.000 0.000
# snapshot number: 2
1.0 1.443747E+00 1.550407E+00 1.514791E+00
2.0 1.403396E+00 2.234155E+00 2.504685E+00
3.0 1.294301E+00 2.264003E+00 2.309210E+00
4.0 2.296760E+00 2.511000E+00 2.881257E-01
5.0 2.297256E+00 2.563023E+00 2.750665E-01
# time, mbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
0.040 5.000 3.000 3.000 3.000 0.000 0.000 0.000
# snapshot number: 3
1.0 1.415620E+00 1.575610E+00 1.522186E+00
2.0 1.444185E+00 2.239136E+00 2.537687E+00
3.0 1.285960E+00 2.272457E+00 2.274797E+00
4.0 2.287877E+00 2.466031E+00 3.080055E-01
5.0 2.301818E+00 2.569355E+00 2.492021E-01
# time, mbeads, dimx, dimy, dimz, shrdx, shrdy, shrdz
0.040 5.000 3.000 3.000 3.000 0.000 0.000 0.000
.. _example:
####################
[Wip] Formatting the HISTORY files of DL_MESO_DPD
####################
.. sidebar:: Software Technical Information
Language
FORTRAN 90
Licence
Documentation Tool
Latex-generated .pdf
Application Documentation
Click to see the manual_
Relevant Training Material
.. Add a link to any relevant training material.
.. contents:: :local:
Purpose of Module
_________________
This module ``format_history.f90`` is a post-processing
utility for DL_MESO_DPD, the Dissipative Particle Dynamics (DPD) code from the DL_MESO_ package.
It converts the trajectory (HISTORY) files from *unformatted* to a
human readable form, (optionally) including explicative comments about all the
quantities. This module is mainly for learning/checking purpouses.
The first aim is to help the user to check that the
system was prepared as intended (e.g., showing all the bead properties and
initial positions, all the bonds etc).
The idea is to use it on small systems when familiarizing with the structure
of input files needed for the simulation.
Secondly, it can be used as a starting point for a user-defined analysis
of trajectories.
.. Mention portability?
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).
.. www.ccp5.ac.uk/DL_MESO
.. M. A. Seaton, R. L. Anderson, S. Metz and W. Smith, Mol. Sim. 39 (10), 796 (2013).
Testing
_______
The present module is compiled with the available Fortran90 compiler, e.g.:
``gfortran -o format.exe format_history.f90``
and the executable must be in the same directory of the HISTORY* files to be
analyzed. To test the module, run the simulation with the *toy* input files given in the following.
(Note that these files contain commented lines as suggestions for further tests.)
For the CONTROL file
.. literalinclude:: ./CONTROL
and for the FIELD file
.. literalinclude:: ./FIELD
After analyzing the trajectories, this output should be printed on the screen
(for both ``lcomm`` and ``lmcheck`` set to ``.TRUE.``):
.. literalinclude:: ./out-co-ch
Source Code
___________
.. literalinclude:: ./format_history.f90
:language: fortran
:linenos:
.. Here are the URL references used
.. _manual: ./manhis.pdf
.. _DL_MESO: http://www.scd.stfc.ac.uk/SCD/support/40694.aspx
.. _ReST: http://docutils.sourceforge.net/docs/user/rst/quickref.html
.. www.ccp5.ac.uk/DL_MESO
../../source/manhis.pdf
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment