1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8 module metaforce 9 10 use precision, only : dp 11 12 implicit none 13 14 integer, save :: nGauss = 0 15 logical, save :: lMetaForce = .false. 16 integer, pointer, save :: nGaussPtr(:,:) => null() 17 real(dp), pointer, save :: GaussK(:) => null() 18 real(dp), pointer, save :: GaussR0(:) => null() 19 real(dp), pointer, save :: GaussZeta(:) => null() 20 21 CONTAINS 22 23 subroutine initmeta() 24C 25C Reads data for metadynamics Gaussians 26C 27C Julian Gale, NRI, Curtin University, Feb. 2006 28C 29 use parallel, only : Node 30 use alloc, only : re_alloc 31 use sys, only : die 32 use fdf 33 34 implicit none 35 36 integer :: k 37 integer :: ni 38 integer :: nr 39 40 type(block_fdf) :: bfdf 41 type(parsed_line), pointer :: pline 42 43C Find if there are any Gaussians specified 44 nGauss = fdf_block_linecount('MetaForce', 'iirrr') 45 lMetaForce = nGauss > 0 46 if ( lMetaForce ) then 47 48C Allocate arrays for Gaussian data 49 call re_alloc( nGaussPtr, 1, 2, 1, nGauss, 'nGaussPtr', 50 & 'initmeta' ) 51 call re_alloc( GaussK, 1, nGauss, 'GaussK', 'initmeta' ) 52 call re_alloc( GaussR0, 1, nGauss, 'GaussR0', 'initmeta' ) 53 call re_alloc( GaussZeta, 1, nGauss, 'GaussZeta', 'initmeta' ) 54 55C Read Gaussians from block 56 lMetaForce = fdf_block('MetaForce', bfdf) 57 58 do while ( fdf_bline(bfdf, pline) ) 59C Read and parse data line 60 ni = fdf_bnintegers(pline) 61 nr = fdf_bnreals(pline) 62 63C Check that correct info is given 64 if ((ni .ge. 2) .and. (nr .ge. 3)) then 65 nGauss = nGauss + 1 66 ! Atom 1 67 nGaussPtr(1,nGauss) = abs(fdf_bintegers(pline,1)) 68 ! Atom 2 69 nGaussPtr(2,nGauss) = abs(fdf_bintegers(pline,2)) 70 ! Gaussian pre-factor for the force 71 ! Unit: Ry 72 GaussK(nGauss) = fdf_breals(pline,1) 73 ! Distance in r - r_0, where r is the distance between 74 ! atom i and j 75 ! Unit: Ang 76 GaussR0(nGauss) = fdf_breals(pline,2)/0.529177_dp 77 ! Gaussian prefactor 78 ! Unit: 1/Bohr**2 79 GaussZeta(nGauss) = fdf_breals(pline,3) 80 endif 81 enddo 82 83 endif 84 85 end subroutine initmeta 86 87 subroutine meta(xa,na,ucell,Emeta,fa,stress,forces,stresses) 88 89 implicit none 90 91C 92C Compute energy and forces due to Gaussians 93C 94C Julian Gale, NRI, Curtin University, Feb. 2006 95C 96 97C Passed variables 98 integer, intent(in) :: na 99 logical, intent(in) :: forces 100 logical, intent(in) :: stresses 101 real(dp), intent(in) :: ucell(3,3) 102 real(dp), intent(in) :: xa(3,na) 103 real(dp), intent(out) :: Emeta 104 real(dp), intent(inout) :: fa(3,na) 105 real(dp), intent(inout) :: stress(3,3) 106 107C Local variables 108 integer :: ii 109 integer :: i 110 integer :: j 111 integer :: k 112 integer :: ixcell(27) 113 integer :: iycell(27) 114 integer :: izcell(27) 115 real(dp) :: dtrm 116 real(dp) :: exptrm 117 real(dp) :: r 118 real(dp) :: r2 119 real(dp) :: r2min 120 real(dp) :: rvol 121 real(dp) :: vol 122 real(dp) :: volcel 123 real(dp) :: x 124 real(dp) :: y 125 real(dp) :: z 126 real(dp) :: xc 127 real(dp) :: yc 128 real(dp) :: zc 129 real(dp) :: xmin 130 real(dp) :: ymin 131 real(dp) :: zmin 132 real(dp) :: xcell(27) 133 real(dp) :: ycell(27) 134 real(dp) :: zcell(27) 135 136 if (stresses) then 137 vol = volcel(ucell) 138 rvol = 1.0_dp/vol 139 endif 140C Initialise energy 141 Emeta = 0.0_dp 142 143C Find cell images 144 call cellimagelist(ucell,xcell,ycell,zcell,ixcell,iycell,izcell) 145 146C Loop over Gaussians 147 do k = 1,nGauss 148 149C Set pointers to atoms 150 i = nGaussPtr(1,k) 151 j = nGaussPtr(2,k) 152 153C Set initial coordinates 154 x = xa(1,j) - xa(1,i) 155 y = xa(2,j) - xa(2,i) 156 z = xa(3,j) - xa(3,i) 157 158C Get coordinates of nearest image 159 r2min = 100000.0_dp 160 do ii = 1,27 161 xc = x + xcell(ii) 162 yc = y + ycell(ii) 163 zc = z + zcell(ii) 164 165C Compute distance 166 r2 = xc*xc + yc*yc + zc*zc 167 if (r2.lt.r2min) then 168 r2min = r2 169 xmin = xc 170 ymin = yc 171 zmin = zc 172 endif 173 174 enddo 175 176 r = sqrt(r2min) 177 178C Compute energy and forces 179 exptrm = exp(-GaussZeta(k)*(r - GaussR0(k))**2) 180 Emeta = Emeta + 0.5_dp*GaussK(k)*exptrm 181 dtrm = GaussK(k)*exptrm*GaussZeta(k)*(r - GaussR0(k))/r 182 if (forces) then 183 fa(1,i) = fa(1,i) - dtrm*xmin 184 fa(2,i) = fa(2,i) - dtrm*ymin 185 fa(3,i) = fa(3,i) - dtrm*zmin 186 fa(1,j) = fa(1,j) + dtrm*xmin 187 fa(2,j) = fa(2,j) + dtrm*ymin 188 fa(3,j) = fa(3,j) + dtrm*zmin 189 endif 190 if (stresses) then 191 dtrm = dtrm*rvol 192 stress(1,1) = stress(1,1) + dtrm*xmin*xmin 193 stress(2,1) = stress(2,1) + dtrm*ymin*xmin 194 stress(3,1) = stress(3,1) + dtrm*zmin*xmin 195 stress(1,2) = stress(1,2) + dtrm*xmin*ymin 196 stress(2,2) = stress(2,2) + dtrm*ymin*ymin 197 stress(3,2) = stress(3,2) + dtrm*zmin*ymin 198 stress(1,3) = stress(1,3) + dtrm*xmin*zmin 199 stress(2,3) = stress(2,3) + dtrm*ymin*zmin 200 stress(3,3) = stress(3,3) + dtrm*zmin*zmin 201 endif 202 203 enddo 204 205 end subroutine meta 206 207 end module metaforce 208