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      subroutine cgvc( na, xa, fa, cell, stress, dxmax,
9     .                 tp, ftol, strtol, varcel, relaxd, usesavecg )
10c ***************************************************************************
11c Variable-cell conjugate-gradient geometry optimization
12c
13c   Energy minimization including atomic coordinates and unit cell vectors.
14c   It allows an external target stress:
15c              %block MD.TargetStress
16c                  3.5  0.0  0.0  0.0  0.0  0.0
17c              %endblock MD.TargetStress
18c   corresponding to xx, yy, zz, xy, xz, yz.
19c   In units of (-MD.TargetPressure)
20c   Default: hydrostatic pressure: -1, -1, -1, 0, 0, 0
21c
22c   Based on E({xa},stress), with {xa} in fractional coor
23c   The gradient of the energy given by {cfa} forces (fractional!)
24c   The fractional coordinates are multiplied by the initial cell vectors
25c   to get them in Bohr for dxmax and preconditioning.
26c
27c Written by E. Artacho. November 1998.
28c ******************************** INPUT ************************************
29c integer na            : Number of atoms in the simulation cell
30c real*8 fa(3,na)       : Atomic forces
31c real*8 dxmax          : Maximum atomic (or lattice vector) displacement
32c real*8 tp             : Target pressure
33c real*8 ftol           : Maximum force tolerance
34c logical varcel        : true if variable cell optimization
35c logical usesavecg     : true if we're using saved CG files.
36c *************************** INPUT AND OUTPUT ******************************
37c real*8 xa(3,na)       : Atomic coordinates
38c                         input: current step; output: next step
39c real*8 cell(3,3)      : Matrix of the vectors defining the unit cell
40c                         input: current step; output: next step
41c                         cell(ix,ja) is the ix-th component of ja-th vector
42c real*8 stress(3,3)    : Stress tensor components
43c real*8 strtol         : Maximum stress tolerance
44c ******************************** OUTPUT ***********************************
45c logical relaxd        : True when converged
46c ***************************************************************************
47
48C
49C  Modules
50C
51      use precision, only : dp
52      use parallel,    only : Node
53      use sys,         only : die
54      use fdf
55
56      use m_mpi_utils, only:  broadcast
57      use units, only: kBar, Ang
58!
59!     Use old-style conjugate-gradient routine
60!     Use the new one only for the Zmatrix case
61!
62      use m_conjgr_old, only: conjgr
63      use alloc,        only: re_alloc, de_alloc
64
65      implicit none
66
67! Subroutine arguments:
68
69      integer, intent(in) :: na
70      real(dp), intent(in) :: fa(3,na), dxmax,
71     .                        tp, ftol, strtol
72      logical, intent(in) :: varcel, usesavecg
73      real(dp), intent(inout) :: xa(3,na), stress(3,3), cell(3,3)
74      logical, intent(out) :: relaxd
75
76c Internal variables and arrays
77
78      real(dp)            :: new_volume, trace, ftol_tmp, volume
79
80      logical           found
81      integer           ia, i, j, n, indi
82
83      real(dp) ::  celli(3,3), sxx, syy, szz, sxy, sxz, syz
84      real(dp) ::  stress_dif(3,3)
85
86      real(dp), pointer :: gxa(:)=>null(), gfa(:)=>null()
87      real(dp), pointer, save :: cgaux(:)=>null()
88
89      type(block_fdf)            :: bfdf
90      type(parsed_line), pointer :: pline=>null()
91
92! Saved internal variables:
93
94      logical, save :: frstme = .true.
95      logical, save :: tarstr = .false.
96      logical, save :: constant_volume
97      real(dp), save :: initial_volume
98
99      real(dp), save :: cgcntr(0:20) = 0.0_dp
100
101      integer, save :: ndeg,
102     .                 linmin
103
104      real(dp), save :: tstres(3,3),
105     .                  modcel(3),
106     .                  precon,
107     .                  strain(3,3),
108     .                  cellin(3,3)
109
110      real(dp) :: volcel
111      external :: volcel
112c ---------------------------------------------------------------------------
113
114      volume = volcel(cell)
115
116C Allocate local memory
117      call re_alloc( gfa, 1, na*3+6, 'gfa', 'cgvc' )
118      call re_alloc( gxa, 1, na*3+6, 'gxa', 'cgvc' )
119      if (.not.associated(cgaux))
120     &  call re_alloc( cgaux, 1, na*6+12, 'cgaux', 'cgvc' )
121
122C If first call to cgvc, check dim and get target stress --------------------
123
124      if ( frstme ) then
125
126        if ( varcel ) then
127
128C Check if we want a constant-volume simulation
129          constant_volume = fdf_get("MD.ConstantVolume", .false.)
130
131C Look for target stress and read it if found, otherwise generate it --------
132
133          tarstr = fdf_block('MD.TargetStress',bfdf)
134
135          if (tarstr) then
136            if (Node.eq.0) then
137              write(6,'(/a,a)') 'cgvc: Reading %block MD.TargetStress',
138     .                          ' (units of MD.TargetPressure).'
139            endif
140            if (.not. fdf_bline(bfdf,pline))
141     .        call die('cgvc: ERROR in MD.TargetStress block')
142            sxx = fdf_bvalues(pline,1)
143            syy = fdf_bvalues(pline,2)
144            szz = fdf_bvalues(pline,3)
145            sxy = fdf_bvalues(pline,4)
146            sxz = fdf_bvalues(pline,5)
147            syz = fdf_bvalues(pline,6)
148            call fdf_bclose(bfdf)
149
150            tstres(1,1) = - sxx * tp
151            tstres(2,2) = - syy * tp
152            tstres(3,3) = - szz * tp
153            tstres(1,2) = - sxy * tp
154            tstres(2,1) = - sxy * tp
155            tstres(1,3) = - sxz * tp
156            tstres(3,1) = - sxz * tp
157            tstres(2,3) = - syz * tp
158            tstres(3,2) = - syz * tp
159          else
160            if (Node.eq.0) then
161              write(6,'(/a,a)') 'cgvc: No target stress found, ',
162     .             'assuming hydrostatic MD.TargetPressure.'
163            endif
164            do i= 1, 3
165              do j= 1, 3
166                tstres(i,j) = 0._dp
167              enddo
168              tstres(i,i) = - tp
169            enddo
170          endif
171
172C Write target stress down --------------------------------------------------
173
174          if (constant_volume) then
175            tstres(:,:) = 0.0_dp
176            if (Node.eq.0) then
177              write(6,"(a)") "***Target stress set to zero " //
178     $            "for constant-volume calculation"
179            endif
180          endif
181          if (Node.eq.0) then
182            write(6,"(/a)") 'cgvc: Target stress (kBar)'
183            write(6,"(a,2x,3f12.3)")
184     .       'cgvc:', tstres(1,1)/kBar, tstres(1,2)/kBar,
185     .       tstres(1,3)/kBar
186            write(6,"(a,2x,3f12.3)")
187     .       'cgvc:', tstres(2,1)/kBar, tstres(2,2)/kBar,
188     .       tstres(2,3)/kBar
189            write(6,"(a,2x,3f12.3)")
190     .       'cgvc:', tstres(3,1)/kBar, tstres(3,2)/kBar,
191     .       tstres(3,3)/kBar
192          endif
193
194
195C Moduli of original cell vectors for fractional coor scaling back to au ---
196
197          do n = 1, 3
198            modcel(n) = 0.0_dp
199            do j = 1, 3
200              modcel(n) = modcel(n) + cell(j,n)*cell(j,n)
201            enddo
202            modcel(n) = dsqrt( modcel(n) )
203          enddo
204
205C Scale factor for strain variables to share magnitude with coordinates
206C ---- (a length in Bohrs typical of bond lengths ..) ------------------
207
208          precon = fdf_get('MD.PreconditionVariableCell',
209     .                     9.4486344d0,'Bohr')
210
211C Dimension of space where E is minimized ------------------------------
212
213          ndeg = na*3 + 6
214
215C Initialize absolute strain and save initial cell vectors -------------
216C Initialization to 1. for numerical reasons, later substracted --------
217
218          strain(1:3,1:3) = 1.0_dp
219          cellin(1:3,1:3) = cell(1:3,1:3)
220          initial_volume = volcel(cellin)
221
222        else
223
224          ndeg = na*3
225
226        endif
227
228C Initialize and read cgaux and cgcntr if present and wanted ---------------
229
230        if (usesavecg) then
231          if (Node.eq.0) then
232            call iocg( 'read', ndeg*2, cgaux, cgcntr, relaxd, found )
233          endif
234          call broadcast(cgaux)
235          call broadcast(cgcntr)
236          call broadcast(relaxd)
237          call broadcast(found)
238
239          if ( found ) then
240            linmin = cgcntr(1)
241! Simple bugfix for fixed cell restarts. In the case of a tightend ftol
242! after a sucessful CG run, relaxed is set to .true. by the read from
243! iocg. This results in the test for convergence being skipped within
244! conjgr. Setting relaxed to .false. here avoids this. In the case of
245! a variable cell relaxed is set below anyway and the converence test
246! is always performed. A better solution (moving the test out into a
247! seperate routine) will follow in a later version.  -- AMW 8/7/2008
248            relaxd = .false.
249          else
250            if (Node.eq.0) then
251              write(6,'(/,a)') 'cgvc: WARNING: CG file not found'
252            endif
253            relaxd = .false.
254            cgcntr(0) = 0
255            linmin = 1
256          endif
257        else
258          relaxd = .false.
259          cgcntr(0) = 0
260          linmin = 1
261        endif
262
263        frstme = .false.
264      endif
265
266C Variable cell -------------------------------------------------------------
267
268      if ( varcel ) then
269
270C Inverse of matrix of cell vectors  (transpose of) ------------------------
271
272        call reclat( cell, celli, 0 )
273
274C Transform coordinates and forces to fractional ----------------------------
275C but scale them again to Bohr by using the (fix) moduli of the original ----
276C lattice vectors (allows using maximum displacement as before) -------------
277C convergence is checked here for input forces as compared with ftol --------
278
279        relaxd = .true.
280        do ia = 1, na
281          do n = 1, 3
282            indi = 3*(ia - 1) + n
283            gxa(indi) = 0.0_dp
284            gfa(indi) = 0.0_dp
285            relaxd = relaxd .and. ( abs(fa(n,ia)) .lt. ftol )
286            do i = 1, 3
287              gxa(indi) = gxa(indi) + celli(i,n) * xa(i,ia) * modcel(n)
288              gfa(indi) = gfa(indi) + cell(i,n) * fa(i,ia) / modcel(n)
289            enddo
290          enddo
291        enddo
292
293C Symmetrizing the stress tensor --------------------------------------------
294
295        do i = 1, 3
296           do j = i+1, 3
297              stress(i,j) = 0.5_dp*( stress(i,j) + stress(j,i) )
298              stress(j,i) = stress(i,j)
299           enddo
300        enddo
301
302C Subtract target stress
303
304        stress_dif = stress - tstres
305!
306!       Take 1/3 of the trace out here if constant-volume needed
307!
308        if (constant_volume) then
309           trace = stress_dif(1,1) + stress_dif(2,2) + stress_dif(3,3)
310           do i=1,3
311              stress_dif(i,i) = stress_dif(i,i) - trace/3.0_dp
312           enddo
313        endif
314
315C Append excess stress and strain to gxa and gfa ------
316C preconditioning: scale stress and strain so as to behave similarly to x,f -
317
318        indi = 3*na
319        do i = 1, 3
320           do j = i, 3
321              indi = indi + 1
322              gfa(indi) = -stress_dif(i,j)*volume/precon
323              gxa(indi) = strain(i,j) * precon
324           enddo
325        enddo
326
327C Check stress convergence --------------------------------------------------
328
329        do i = 1, 3
330           do j = 1, 3
331              relaxd = relaxd .and.
332     .          ( abs(stress_dif(i,j)) .lt. abs(strtol) )
333           enddo
334        enddo
335
336C Call conjugate gradient minimization --------------------------------------
337
338        if ( .not. relaxd ) then
339           ftol_tmp = 0.0_dp
340!!           do i=1,ndeg
341!!              print *, "gxa, gfa ", i, gxa(i), gfa(i)
342!!           enddo
343           call conjgr(ndeg,gxa,gfa,dxmax, ftol_tmp ,cgcntr,cgaux )
344        endif
345
346
347C Fixed cell ----------------------------------------------------------------
348
349      else
350
351         if (.not. relaxd) then
352!!            do i=1,ndeg
353!!               print *, "gxa, gfa ", i, gxa(i), gfa(i)
354!!            enddo
355            call conjgr( 3*na, xa, fa, dxmax, ftol, cgcntr, cgaux )
356            relaxd = (int(cgcntr(0)) .eq. 0)           !! ???
357         endif
358
359      endif
360
361C Checking line minimizations and convergence -------------------------------
362
363      if (nint(cgcntr(1)) .ne. linmin) then
364        if (Node.eq.0) then
365          write(6,'(/a,i4,a,f10.4)')
366     .      'cgvc: Finished line minimization ', linmin,
367     .      '.  Mean atomic displacement =', cgcntr(18)/sqrt(dble(na))
368        endif
369        linmin = nint(cgcntr(1))
370      endif
371
372C Transform back if variable cell -------------------------------------------
373
374      if ( varcel .and. (.not. relaxd) ) then
375
376C New cell ------------------------------------------------------------------
377
378        indi = 3*na
379        do i = 1, 3
380           do j = i, 3
381              indi = indi + 1
382              strain(i,j) = gxa(indi) / precon
383              strain(j,i) = strain(i,j)
384           enddo
385        enddo
386
387        cell = cellin + matmul(strain-1.0_dp,cellin)
388        if (constant_volume) then
389           new_volume = volcel(cell)
390           if (Node.eq.0) write(6,"(a,f12.4)")
391     $          "Volume before coercion: ",  new_volume/Ang**3
392           cell =  cell * (initial_volume/new_volume)**(1.0_dp/3.0_dp)
393        endif
394
395C Output fractional coordinates to cartesian Bohr, and copy to xa -----------
396
397        do ia = 1, na
398          do i = 1, 3
399            xa(i,ia) = 0.0_dp
400            do n = 1, 3
401              indi = 3*(ia - 1) + n
402              xa(i,ia) = xa(i,ia) + cell(i,n) * gxa(indi) / modcel(n)
403            enddo
404          enddo
405        enddo
406
407      endif
408
409C Save cgaux ----------------------------------------------------------------
410
411      if (Node.eq.0)
412     .  call iocg( 'write', ndeg*2, cgaux, cgcntr, relaxd, found )
413
414C Deallocate local memory
415      call de_alloc( gxa, 'gxa', 'cgvc' )
416      call de_alloc( gfa, 'gfa', 'cgvc' )
417
418      return
419      end
420
421