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 m_cell_fire_optim
9!
10!     Fire geometry optimization
11!
12      use precision, only : dp
13      use sys, only: die
14      use alloc
15      use fdf, only: fdf_get
16
17      use m_fire
18      use siesta_options, only: dt
19
20      use m_memory, only: memory, mem_stat
21      use parallel, only : ionode
22
23      implicit none
24
25      public :: cell_fire_optimizer
26      private
27
28      CONTAINS
29      subroutine cell_fire_optimizer( na, xa, cell, stress,
30     $           tp, strtol, varcel, relaxd)
31c ***************************************************************************
32c Fire geometry optimization (Cell Only)
33c
34c Written by Alberto Garcia, April 2007
35
36c ******************************** INPUT ************************************
37c na                    : number of atoms (will need to keep fractional coords)
38c real*8 tp             : Target pressure
39c logical varcel        : true if variable cell optimization
40c *************************** INPUT AND OUTPUT ******************************
41c real*8 cell(3,3)      : Matrix of the vectors defining the unit cell
42c                         input: current step; output: next step
43c                         cell(ix,ja) is the ix-th component of ja-th vector
44c real*8 xa(3,na)       : atomic coordinates
45c real*8 stress(3,3)    : Stress tensor components
46c real*8 strtol         : Maximum stress tolerance
47c ******************************** OUTPUT ***********************************
48c logical relaxd        : True when converged
49c ***************************************************************************
50
51C
52C  Modules
53C
54
55      implicit none
56! Subroutine arguments:
57
58      real(dp), intent(in) :: tp,  strtol
59      logical, intent(in) :: varcel
60      real(dp), intent(inout) :: stress(3,3), cell(3,3)
61      integer, intent(in)  :: na
62      real(dp), intent(inout) :: xa(3,na)
63      logical, intent(out) :: relaxd
64
65c Internal variables and arrays
66
67      real(dp)       :: volume, new_volume, trace
68
69      integer           i, j, indi
70
71      real(dp) ::  celli(3,3)
72      real(dp) ::  stress_dif(3,3)
73
74      real(dp), dimension(:), allocatable       :: gxa, gfa
75      real(dp), dimension(:), pointer       :: deltamax
76
77
78! Saved internal variables:
79
80      logical, save :: frstme = .true.
81      logical, save :: constant_volume
82      real(dp), save :: initial_volume
83
84
85      real(dp), save :: target_stress(3,3),
86     .                  precon,
87     .                  strain(3,3),
88     .                  cellin(3,3)
89
90      type(fire_t), save  :: b
91      integer, save  :: numel
92
93      logical, save  :: fire_debug
94      real(dp), save :: fire_mass
95      real(dp)       :: fire_dt, fire_dt_inc,
96     $                  fire_dt_dec, fire_alphamax,
97     $                  fire_alpha_factor, fire_dtmax
98      integer        :: fire_nmin
99      real(dp), parameter ::  dstrain_max = 0.1_dp
100
101      real(dp) :: volcel
102      external :: volcel
103c ---------------------------------------------------------------------------
104
105      volume = volcel(cell)
106
107      if ( frstme ) then
108         fire_mass = fdf_get("MD.FIRE.Mass", 1.0_dp)
109         fire_dt = fdf_get("MD.FIRE.TimeStep", dt)
110         fire_dt_inc = fdf_get("MD.FIRE.TimeInc", FIRE_DEF_dt_inc)
111         fire_dt_dec = fdf_get("MD.FIRE.TimeDec", FIRE_DEF_dt_dec)
112         fire_nmin = fdf_get("MD.FIRE.Nmin", FIRE_DEF_nmin)
113         fire_alphamax = fdf_get("MD.FIRE.AlphaMax", FIRE_DEF_alphamax)
114         fire_alpha_factor = fdf_get("MD.FIRE.AlphaFactor",
115     &        FIRE_DEF_alpha_factor)
116         fire_dtmax = fdf_get("MD.FIRE.MaxTimeStep", FIRE_DEF_dtmax)
117         fire_dt = fdf_get("MD.FIRE.TimeStep", dt)
118         fire_debug = fdf_get("MD.FIRE.Debug", .false.)
119
120
121         if (varcel ) then
122            numel =  6
123         else
124            call die("no varcel?")
125         endif
126         if (Ionode) then
127           write(6,'(a,i6)') "Cell_Fire_optim: No of elements: ",
128     $                numel
129         endif
130
131         call fire_setup(b, n=numel, dt=fire_dt,
132     $                   debug=fire_debug,
133     $                   dt_inc=fire_dt_inc, dt_dec=fire_dt_dec,
134     $                   alphamax=fire_alphamax,
135     $                   alpha_factor=fire_alpha_factor,
136     $                   nmin=fire_nmin)
137
138        if ( varcel ) then
139
140C Check if we want a constant-volume simulation
141          constant_volume = fdf_get("MD.ConstantVolume", .false.)
142
143          call get_target_stress(tp,target_stress)
144          if (constant_volume) target_stress = 0.0_dp
145
146
147C Scale factor for strain variables to share magnitude with coordinates
148C ---- (a length in Bohrs typical of bond lengths ..) ------------------
149
150          precon = fdf_get('MD.PreconditionVariableCell',
151     .                           9.4486344d0,'Bohr')
152
153C Initialize absolute strain and save initial cell vectors -------------
154C Initialization to 1. for numerical reasons, later substracted --------
155
156          strain(1:3,1:3) = 1.0_dp
157          cellin(1:3,1:3) = cell(1:3,1:3)
158          initial_volume = volcel(cellin)
159
160        endif
161
162        frstme = .false.
163      endif                     ! First call
164
165C Variable cell -------------------------------------------------------------
166
167      if ( varcel ) then
168
169         allocate(gfa(numel), stat=mem_stat)
170         call memory('A','D',numel,'Fire_optim',
171     $        stat=mem_stat,id="gfa")
172         allocate(gxa(numel), stat=mem_stat)
173         call memory('A','D',numel,'Fire_optim',
174     $        stat=mem_stat,id="gxa")
175        nullify( deltamax )
176        call re_alloc( deltamax, 1, numel, name='deltamax',
177     &                 routine='fire_optimizer' )
178
179
180        relaxd = .true.
181
182C Symmetrizing the stress tensor
183
184        do i = 1, 3
185           do j = i+1, 3
186              stress(i,j) = 0.5_dp*( stress(i,j) + stress(j,i) )
187              stress(j,i) = stress(i,j)
188           enddo
189        enddo
190
191C Subtract target stress
192
193        stress_dif = stress - target_stress
194!
195!       Take 1/3 of the trace out here if constant-volume needed
196!
197        if (constant_volume) then
198           trace = stress_dif(1,1) + stress_dif(2,2) + stress_dif(3,3)
199           do i=1,3
200              stress_dif(i,i) = stress_dif(i,i) - trace/3.0_dp
201           enddo
202        endif
203
204C Append excess stress and strain to gxa and gfa ------
205C preconditioning: scale stress and strain so as to behave similarly to x,f -
206
207        indi = 0
208        do i = 1, 3
209           do j = i, 3
210              indi = indi + 1
211              gfa(indi) = -stress_dif(i,j)*volume/precon
212              gxa(indi) = strain(i,j) * precon
213              deltamax(indi) = dstrain_max
214           enddo
215        enddo
216
217C Check stress convergence --------------------------------------------------
218
219        do i = 1, 3
220           do j = 1, 3
221              relaxd = relaxd .and.
222     .          ( abs(stress_dif(i,j)) .lt. abs(strtol) )
223           enddo
224        enddo
225
226        if (relaxd) RETURN
227
228C Call Fire step
229
230        call fire_step(b,gfa,gxa,deltamax)
231
232      endif
233
234C Transform back if variable cell
235
236      if ( varcel ) then
237
238      ! New cell
239
240        indi = 0
241        do i = 1, 3
242           do j = i, 3
243              indi = indi + 1
244              strain(i,j) = gxa(indi) / precon
245              strain(j,i) = strain(i,j)
246           enddo
247        enddo
248
249        ! Inverse of matrix of cell vectors  (transpose of)
250        call reclat( cell, celli, 0 )
251
252!       Update cell
253
254        cell = cellin + matmul(strain-1.0_dp,cellin)
255        if (constant_volume) then
256           new_volume = volcel(cell)
257           cell =  cell * (initial_volume/new_volume)**(1.0_dp/3.0_dp)
258        endif
259
260        call rescale_coordinates(na,xa, celli_oldcell=celli,
261     $                                  new_cell=cell)
262
263
264      ! Deallocate local memory
265
266        deallocate (gxa, stat=mem_stat)
267        call memory('D','D',numel,'Fire_optim',
268     $       stat=mem_stat,id="gxa")
269        deallocate (gfa, stat=mem_stat)
270        call memory('D','D',numel,'Fire_optim',
271     $       stat=mem_stat,id="gfa")
272        call de_alloc( deltamax, name='deltamax' )
273
274      endif ! variable cell
275
276      end subroutine cell_fire_optimizer
277!
278!--------------------------------------------------------------
279
280      subroutine rescale_coordinates(na,xa,
281     $                            celli_oldcell,new_cell)
282      use precision, only : dp
283      use zmatrix
284
285      integer, intent(in)     :: na
286      real(dp), intent(inout) :: xa(3,na)
287      real(dp), intent(in)    :: celli_oldcell(3,3)
288      real(dp), intent(in)    :: new_cell(3,3)
289
290
291      real(dp), dimension(3) :: r, frac
292      integer  :: ifirst, imol, icart, i, j
293
294      !  NOTE: We have to be careful here if using a Zmatrix
295      if (lUseZmatrix) then
296
297        !     re-scale only the positions of the leading atoms
298        !     in each molecule,
299        !     plus any cartesian block,
300        !     and recompute the cartesian coordinate array
301        !
302           do imol = 1, nZmol
303              ifirst = nZmolStartAtom(imol)
304              r(1:3) = Zmat(3*(ifirst-1)+1:3*(ifirst-1)+3)
305              frac(1:3) = matmul(transpose(celli_oldcell),r(1:3))
306              r(1:3) = matmul(new_cell,frac(1:3))
307              Zmat(3*(ifirst-1)+1:3*(ifirst-1)+3) = r(1:3)
308           enddo
309           do icart = 1, nZcart
310             ! Process cartesian block
311              ifirst = nZcartStartAtom(icart)
312              do j = ifirst, ifirst + nZcartAtoms(icart) - 1
313                 r(1:3) = Zmat(3*(j-1)+1:3*(j-1)+3)
314                 frac(1:3) = matmul(transpose(celli_oldcell),r(1:3))
315                 Zmat(3*(j-1)+1:3*(j-1)+3) = matmul(new_cell,frac(1:3))
316              enddo
317           enddo
318           call Zmat_to_Cartesian(xa)
319        else
320           ! No Zmatrix
321           ! Rescale coordinates for all atoms
322           do i = 1, na
323              xa(:,i) = matmul(transpose(celli_oldcell),xa(:,i))
324              xa(:,i) = matmul(new_cell,xa(:,i))
325           enddo
326
327        endif
328
329      end subroutine rescale_coordinates
330
331
332      end module m_cell_fire_optim
333
334