forked from lijiext/lammps
216 lines
9.0 KiB
Fortran
216 lines
9.0 KiB
Fortran
! This performs the calculations necessary to run a Lennard-Jones (LJ)
|
|
! simulation.
|
|
!
|
|
! Copyright (C) 2013, Joshua More and Michele Ceriotti
|
|
!
|
|
! Permission is hereby granted, free of charge, to any person obtaining
|
|
! a copy of this software and associated documentation files (the
|
|
! "Software"), to deal in the Software without restriction, including
|
|
! without limitation the rights to use, copy, modify, merge, publish,
|
|
! distribute, sublicense, and/or sell copies of the Software, and to
|
|
! permit persons to whom the Software is furnished to do so, subject to
|
|
! the following conditions:
|
|
!
|
|
! The above copyright notice and this permission notice shall be included
|
|
! in all copies or substantial portions of the Software.
|
|
!
|
|
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
|
! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
|
|
! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
|
|
! CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
|
|
! TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
|
|
! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
|
|
!
|
|
!
|
|
! This contains the functions that calculate the potential, forces and
|
|
! virial tensor of a single-component LJ system.
|
|
! Includes functions which calculate the long-range correction terms for a
|
|
! simulation with a sharp nearest-neighbour cut-off.
|
|
!
|
|
! Functions:
|
|
! LJ_functions: Calculates the LJ pair potential and the magnitude of the
|
|
! forces acting on a pair of atoms.
|
|
! LJ_fij: Calculates the LJ pair potential and force vector for the
|
|
! interaction of a pair of atoms.
|
|
! LJ_longrange: Calculates the long range correction to the potential
|
|
! and virial.
|
|
! LJ_getall: Calculates the potential and virial of the system and the
|
|
! forces acting on all the atoms.
|
|
|
|
MODULE LJ
|
|
USE DISTANCE
|
|
IMPLICIT NONE
|
|
|
|
DOUBLE PRECISION, PARAMETER :: four_tau_by_3 = 8.3775804095727811d0
|
|
|
|
CONTAINS
|
|
|
|
SUBROUTINE LJ_functions(sigma, eps, r, pot, force)
|
|
! Calculates the magnitude of the LJ force and potential between
|
|
! a pair of atoms at a given distance from each other.
|
|
!
|
|
! Args:
|
|
! sigma: The LJ distance parameter.
|
|
! eps: The LJ energy parameter.
|
|
! r: The separation of the atoms.
|
|
! pot: The LJ interaction potential.
|
|
! force: The magnitude of the LJ force.
|
|
|
|
DOUBLE PRECISION, INTENT(IN) :: sigma
|
|
DOUBLE PRECISION, INTENT(IN) :: eps
|
|
DOUBLE PRECISION, INTENT(IN) :: r
|
|
DOUBLE PRECISION, INTENT(OUT) :: pot
|
|
DOUBLE PRECISION, INTENT(OUT) :: force
|
|
|
|
DOUBLE PRECISION sigma_by_r6
|
|
|
|
sigma_by_r6 = sigma/r
|
|
sigma_by_r6 = sigma_by_r6*sigma_by_r6*sigma_by_r6
|
|
sigma_by_r6 = sigma_by_r6*sigma_by_r6
|
|
|
|
pot = 4*eps*(sigma_by_r6*(sigma_by_r6 - 1))
|
|
force = 24*eps*(sigma_by_r6*(2*sigma_by_r6 - 1)/r)
|
|
|
|
END SUBROUTINE
|
|
|
|
SUBROUTINE LJ_fij(sigma, eps, rij, r, pot, fij)
|
|
! This calculates the LJ potential energy and the magnitude and
|
|
! direction of the force acting on a pair of atoms.
|
|
!
|
|
! Args:
|
|
! sigma: The LJ distance parameter.
|
|
! eps: The LJ energy parameter.
|
|
! rij: The vector joining the two atoms.
|
|
! r: The separation of the two atoms.
|
|
! pot: The LJ interaction potential.
|
|
! fij: The LJ force vector.
|
|
|
|
DOUBLE PRECISION, INTENT(IN) :: sigma
|
|
DOUBLE PRECISION, INTENT(IN) :: eps
|
|
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: rij
|
|
DOUBLE PRECISION, INTENT(IN) :: r
|
|
DOUBLE PRECISION, INTENT(OUT) :: pot
|
|
DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: fij
|
|
|
|
DOUBLE PRECISION f_tot
|
|
|
|
CALL LJ_functions(sigma, eps, r, pot, f_tot)
|
|
fij = f_tot*rij/r
|
|
|
|
END SUBROUTINE
|
|
|
|
SUBROUTINE LJ_longrange(rc, sigma, eps, natoms, volume, pot_lr, vir_lr)
|
|
! Calculates the long range correction to the total potential and
|
|
! virial pressure.
|
|
!
|
|
! Uses the tail correction for a sharp cut-off, with no smoothing
|
|
! function, as derived in Martyna and Hughes, Journal of Chemical
|
|
! Physics, 110, 3275, (1999).
|
|
!
|
|
! Args:
|
|
! rc: The cut-off radius.
|
|
! sigma: The LJ distance parameter.
|
|
! eps: The LJ energy parameter.
|
|
! natoms: The number of atoms in the system.
|
|
! volume: The volume of the system box.
|
|
! pot_lr: The tail correction to the LJ interaction potential.
|
|
! vir_lr: The tail correction to the LJ virial pressure.
|
|
|
|
DOUBLE PRECISION, INTENT(IN) :: rc
|
|
DOUBLE PRECISION, INTENT(IN) :: sigma
|
|
DOUBLE PRECISION, INTENT(IN) :: eps
|
|
INTEGER, INTENT(IN) :: natoms
|
|
DOUBLE PRECISION, INTENT(IN) :: volume
|
|
DOUBLE PRECISION, INTENT(OUT) :: pot_lr
|
|
DOUBLE PRECISION, INTENT(OUT) :: vir_lr
|
|
|
|
DOUBLE PRECISION sbyr, s3byr3, s6byr3, s6byr6, prefactor
|
|
|
|
sbyr = sigma/rc
|
|
s3byr3 = sbyr*sbyr*sbyr
|
|
s6byr6 = s3byr3*s3byr3
|
|
prefactor = four_tau_by_3*natoms*natoms*eps/volume
|
|
prefactor = prefactor*s3byr3*sigma*sigma*sigma
|
|
|
|
pot_lr = prefactor*(s6byr6/3-1)
|
|
vir_lr = prefactor*(s6byr6-1) + pot_lr
|
|
|
|
END SUBROUTINE
|
|
|
|
SUBROUTINE LJ_getall(rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list, n_list, pot, forces, virial)
|
|
! Calculates the LJ potential energy and virial and the forces
|
|
! acting on all the atoms.
|
|
!
|
|
! Args:
|
|
! rc: The cut-off radius.
|
|
! sigma: The LJ distance parameter.
|
|
! eps: The LJ energy parameter.
|
|
! natoms: The number of atoms in the system.
|
|
! atoms: A vector holding all the atom positions.
|
|
! cell_h: The simulation box cell vector matrix.
|
|
! cell_ih: The inverse of the simulation box cell vector matrix.
|
|
! index_list: A array giving the last index of n_list that
|
|
! gives the neighbours of a given atom.
|
|
! n_list: An array giving the indices of the atoms that neighbour
|
|
! the atom determined by index_list.
|
|
! pot: The total potential energy of the system.
|
|
! forces: An array giving the forces acting on all the atoms.
|
|
! virial: The virial tensor, not divided by the volume.
|
|
|
|
DOUBLE PRECISION, INTENT(IN) :: rc
|
|
DOUBLE PRECISION, INTENT(IN) :: sigma
|
|
DOUBLE PRECISION, INTENT(IN) :: eps
|
|
INTEGER, INTENT(IN) :: natoms
|
|
DOUBLE PRECISION, DIMENSION(natoms,3), INTENT(IN) :: atoms
|
|
DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_h
|
|
DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_ih
|
|
INTEGER, DIMENSION(natoms), INTENT(IN) :: index_list
|
|
INTEGER, DIMENSION(natoms*(natoms-1)/2), INTENT(IN) :: n_list
|
|
DOUBLE PRECISION, INTENT(OUT) :: pot
|
|
DOUBLE PRECISION, DIMENSION(natoms,3), INTENT(OUT) :: forces
|
|
DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: virial
|
|
|
|
INTEGER i, j, k, l, start
|
|
DOUBLE PRECISION, DIMENSION(3) :: fij, rij
|
|
DOUBLE PRECISION r2, pot_ij, pot_lr, vir_lr, volume
|
|
|
|
forces = 0.0d0
|
|
pot = 0.0d0
|
|
virial = 0.0d0
|
|
|
|
start = 1
|
|
|
|
DO i = 1, natoms - 1
|
|
! Only loops over the neigbour list, not all the atoms.
|
|
DO j = start, index_list(i)
|
|
CALL vector_separation(cell_h, cell_ih, atoms(i,:), atoms(n_list(j),:), rij, r2)
|
|
IF (r2 < rc*rc) THEN ! Only calculates contributions between neighbouring particles.
|
|
CALL LJ_fij(sigma, eps, rij, sqrt(r2), pot_ij, fij)
|
|
|
|
forces(i,:) = forces(i,:) + fij
|
|
forces(n_list(j),:) = forces(n_list(j),:) - fij
|
|
pot = pot + pot_ij
|
|
DO k = 1, 3
|
|
DO l = k, 3
|
|
! Only the upper triangular elements calculated.
|
|
virial(k,l) = virial(k,l) + fij(k)*rij(l)
|
|
ENDDO
|
|
ENDDO
|
|
ENDIF
|
|
ENDDO
|
|
start = index_list(i) + 1
|
|
ENDDO
|
|
|
|
! Assuming an upper-triangular vector matrix for the simulation box.
|
|
volume = cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
|
|
CALL LJ_longrange(rc, sigma, eps, natoms, volume, pot_lr, vir_lr)
|
|
pot = pot + pot_lr
|
|
DO k = 1, 3
|
|
virial(k,k) = virial(k,k) + vir_lr
|
|
ENDDO
|
|
|
|
END SUBROUTINE
|
|
|
|
END MODULE
|