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