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_zm_broyden_optim
9      public :: zm_broyden_optimizer
10      private
11
12      CONTAINS
13!
14      subroutine zm_broyden_optimizer( na, xa, cell, stress, tp, strtol,
15     &                                 varcel, relaxd )
16c ******************************** INPUT ************************************
17c integer na            : Number of atoms in the simulation cell
18c real*8 stress(3,3)    : Stress tensor components
19c real*8 tp             : Target pressure
20c real*8 strtol         : Maximum stress tolerance
21c logical varcel        : true if variable cell optimization
22c *************************** INPUT AND OUTPUT ******************************
23c real*8 xa(3,na)       : Atomic coordinates
24c                         input: current step; output: next step
25c real*8 cell(3,3)      : Matrix of the vectors defining the unit cell
26c                         input: current step; output: next step
27c                         cell(ix,ja) is the ix-th component of ja-th vector
28c real*8 stress(3,3)    : Stress tensor components
29c real*8 strtol         : Maximum stress tolerance
30c ******************************** OUTPUT ***********************************
31c logical relaxd        : True when converged
32c ***************************************************************************
33
34C
35C  Modules
36C
37      use precision,   only : dp
38      use parallel,    only : Node, IONode
39      use sys,         only : die
40      use fdf
41      use units, only: kBar, Ang
42      use m_broyddj_nocomm, only: broyden_t
43      use m_broyddj_nocomm, only: broyden_reset, broyden_step,
44     $                     broyden_destroy, broyden_init,
45     $                     broyden_is_setup
46
47      use zmatrix,     only : VaryZmat, Zmat
48      use zmatrix,     only : CartesianForce_to_ZmatForce
49      use zmatrix,     only : ZmatForce,ZmatForceVar
50      use zmatrix,     only : iZmattoVars,ZmatType
51      use zmatrix,     only : Zmat_to_Cartesian
52      use zmatrix,     only : coeffA, coeffB, iNextDept
53      use Zmatrix,     only : ZmatForceTolLen, ZmatForceTolAng
54      use Zmatrix,     only : ZmatMaxDisplLen, ZmatMaxDisplAng
55
56      implicit none
57
58C Subroutine arguments:
59      integer,  intent(in)   :: na
60      logical,  intent(in)   :: varcel
61      logical,  intent(out)  :: relaxd
62      real(dp), intent(in)   :: tp
63      real(dp), intent(inout):: xa(3,na), stress(3,3), strtol, cell(3,3)
64
65C Internal variables and arrays
66      integer             :: ia, i, j, n, indi,indi1,vi,k
67      real(dp)            :: volume, new_volume, trace
68      real(dp)            :: sxx, syy, szz, sxy, sxz, syz
69      real(dp)            :: celli(3,3)
70      real(dp)            :: stress_dif(3,3)
71      real(dp), allocatable :: gxa(:), gfa(:)
72      real(dp)            :: force, force1
73
74C Saved internal variables:
75      integer,       save :: ndeg
76      logical,       save :: frstme = .true.
77      logical,       save :: tarstr = .false.
78      logical,       save :: constant_volume
79      real(dp),      save :: initial_volume
80
81      real(dp),      save :: tstres(3,3)
82      real(dp),      save :: cellin(3,3)
83      real(dp),      save :: modcel(3)
84      real(dp),      save :: precon
85      real(dp),      save :: strain(3,3)  ! Special treatment !!
86
87      real(dp), allocatable :: ftoln(:), dxmaxn(:), rnew(:)
88
89      type(broyden_t), save  :: br
90      logical, save           :: initialization_done = .false.
91
92      type(block_fdf)            :: bfdf
93      type(parsed_line), pointer :: pline
94
95      real(dp), save :: jinv0
96      integer, save  :: maxit
97      logical, save  :: cycle_on_maxit, variable_weight
98      logical, save  :: broyden_debug
99      real(dp)       :: weight
100
101      integer        :: ndegi, ndi
102
103      real(dp), external :: volcel
104
105c ---------------------------------------------------------------------------
106
107      volume = volcel(cell)
108
109      if (.not. initialization_done) then
110
111        maxit           = fdf_get("MD.Broyden.History.Steps",5)
112        cycle_on_maxit  = fdf_get("MD.Broyden.Cycle.On.Maxit",.true.)
113        variable_weight = fdf_get("MD.Broyden.Variable.Weight",.false.)
114        broyden_debug   = fdf_get("MD.Broyden.Debug",.false.)
115        jinv0           = fdf_get("MD.Broyden.Initial.Inverse.Jacobian",
116     $                            1.0_dp)
117
118        if (ionode) then
119          write(6,*)
120          write(6,'(a,i3)')
121     $         "Broyden_optim: max_history for broyden: ", maxit
122          write(6,'(a,l1)')
123     $         "Broyden_optim: cycle on maxit: ", cycle_on_maxit
124          if (variable_weight) write(6,'(a)')
125     $         "Broyden_optim: Variable weight not implemented yet"
126          write(6,'(a,f8.4)')
127     $         "Broyden_optim: initial inverse jacobian: ", jinv0
128          write(6,*)
129        endif
130
131        call broyden_init(br,debug=broyden_debug)
132        initialization_done = .true.
133
134      endif
135
136      if ( frstme ) then
137
138        ! Find number of variables
139        ndeg = 0
140        do ia = 1,na
141           do n = 1,3
142              indi = 3*(ia-1) + n
143              if (VaryZmat(indi)) then
144                 ndeg = ndeg + 1
145              endif
146           enddo
147        enddo
148        if ( varcel ) then
149           ndeg = ndeg + 6      ! Including stress
150        endif
151        if (Ionode) then
152           write(6,'(a,i6)') "Broyden_optim: No of elements: ", ndeg
153        endif
154
155        if ( varcel ) then
156
157          constant_volume = fdf_get("MD.ConstantVolume", .false.)
158
159          tarstr = fdf_block('MD.TargetStress',bfdf)
160
161          ! Look for target stress and read it if found,
162          ! otherwise generate it
163          if (tarstr) then
164            if (IOnode) then
165              write(6,'(/a,a)')
166     $          'zm_broyden_optimizer: Reading %block MD.TargetStress',
167     .          ' (units of MD.TargetPressure).'
168            endif
169            if (.not. fdf_bline(bfdf,pline))
170     $        call die('zm_broyden_optimizer: ERROR in ' //
171     .                 'MD.TargetStress block')
172            sxx = fdf_bvalues(pline,1)
173            syy = fdf_bvalues(pline,2)
174            szz = fdf_bvalues(pline,3)
175            sxy = fdf_bvalues(pline,4)
176            sxz = fdf_bvalues(pline,5)
177            syz = fdf_bvalues(pline,6)
178            call fdf_bclose(bfdf)
179
180            tstres(1,1) = - sxx * tp
181            tstres(2,2) = - syy * tp
182            tstres(3,3) = - szz * tp
183            tstres(1,2) = - sxy * tp
184            tstres(2,1) = - sxy * tp
185            tstres(1,3) = - sxz * tp
186            tstres(3,1) = - sxz * tp
187            tstres(2,3) = - syz * tp
188            tstres(3,2) = - syz * tp
189          else
190            write(6,'(/a,a)')
191     $        'zm_broyden_optimizer: No target stress found, ',
192     .        'assuming hydrostatic MD.TargetPressure.'
193            do i = 1, 3
194              do j = 1, 3
195                tstres(i,j) = 0.0_dp
196              enddo
197              tstres(i,i) = - tp
198            enddo
199          endif
200
201          if (constant_volume) then
202            tstres(:,:) = 0.0_dp
203            if (IOnode) then
204              write(6,"(a)") "***Target stress set to zero " //
205     $                       "for constant-volume calculation"
206            endif
207          endif
208          if (IOnode) then
209            write(6,"(/a)") 'cgvc_zmatrix: Target stress (kBar)'
210            do i=1, 3
211              write(6,"(a,2x,3f12.3)")
212     .          'cgvc_zmatrix:', (tstres(i,j)/kBar,j=1,3)
213            enddo
214          endif
215
216C Moduli of original cell vectors for fractional coor scaling back to au ---
217          do n = 1, 3
218            modcel(n) = 0.0_dp
219            do j = 1, 3
220              modcel(n) = modcel(n) + cell(j,n)*cell(j,n)
221            enddo
222            modcel(n) = dsqrt( modcel(n) )
223          enddo
224
225C Scale factor for strain variables to share magnitude with coordinates
226C ---- (a length in Bohrs typical of bond lengths ..)
227
228          ! AG: This could better be volume^(1/3) by default
229          precon = fdf_get('MD.PreconditionVariableCell',
230     .                     9.4486344d0,'Bohr')
231
232C Initialize absolute strain and save initial cell vectors
233C Initialization to 1 for numerical reasons, later substracted
234
235          strain(1:3,1:3) = 1.0_dp
236          cellin(1:3,1:3) = cell(1:3,1:3)
237          initial_volume = volcel(cellin)
238
239        endif     ! varcel
240
241        frstme = .false.
242      endif                 ! First time
243
244C Allocate local memory
245      allocate(gfa(1:ndeg),gxa(1:ndeg))
246      allocate(ftoln(1:ndeg),dxmaxn(1:ndeg))
247
248!      print *, "zm_broyden: cell in : ", cell
249!      print *, "zm_broyden: strain in : ", strain
250
251      if ( varcel ) then
252
253        ! Inverse of matrix of cell vectors  (transpose of)
254        call reclat( cell, celli, 0 )
255
256C Transform coordinates and forces to fractional ----------------------------
257C but scale them again to Bohr by using the (fix) moduli of the original ----
258C lattice vectors (allows using maximum displacement as before) -------------
259C convergence is checked here for input forces and compared with ftoln ------
260
261        relaxd = .true.
262        ndegi = 0
263        ! Loop over degrees of freedom, scaling
264        ! only cartesian coordinates to fractional
265          do ia = 1,na
266            do n = 1,3
267              indi = 3*(ia-1) + n
268              if (VaryZmat(indi)) then
269                ndegi = ndegi + 1
270                if (ZmatType(indi).gt.1) then
271                  ftoln(ndegi) = ZmatForceTolLen
272                  dxmaxn(ndegi) = ZmatMaxDisplLen
273                else
274                  ftoln(ndegi) = ZmatForceTolAng
275                  dxmaxn(ndegi) = ZmatMaxDisplAng
276                endif
277                vi = iZmattoVars(indi)
278                if (vi.eq.0) then
279                  force = ZmatForce(indi)
280                else
281                  force = ZmatForceVar(vi)
282                endif
283                relaxd=relaxd.and.(dabs(force).lt.ftoln(ndegi))
284                if (ZmatType(indi).gt.2) then
285C Cartesian coordinate
286                  gxa(ndegi) = 0.0_dp
287                  gfa(ndegi) = 0.0_dp
288                  do i = 1,3
289                    indi1 = 3*(ia-1)+i
290                    gxa(ndegi) = gxa(ndegi)+
291     .                          celli(i,n)*Zmat(indi1)*modcel(n)
292                    if (i.eq.n) then
293                      force1 = force
294                    else
295                      force1 = ZmatForce(indi1)
296                    endif
297                    gfa(ndegi) = gfa(ndegi)+
298     .                          cell(i,n)*force1/modcel(n)
299                  enddo
300                else
301                  gxa(ndegi) = Zmat(indi)
302                  gfa(ndegi) = force
303                endif
304              endif
305            enddo
306          enddo
307
308C Symmetrizing the stress tensor --------------------------------------------
309        do i = 1,3
310          do j = i+1,3
311            stress(i,j) = 0.5_dp*( stress(i,j) + stress(j,i) )
312            stress(j,i) = stress(i,j)
313          enddo
314        enddo
315
316C Subtract target stress
317
318        stress_dif = stress - tstres
319!
320!       Take 1/3 of the trace out here if constant-volume needed
321!
322        if (constant_volume) then
323           trace = stress_dif(1,1) + stress_dif(2,2) + stress_dif(3,3)
324           do i=1,3
325              stress_dif(i,i) = stress_dif(i,i) - trace/3.0_dp
326           enddo
327        endif
328
329C Append stress (substracting target stress) and strain to gxa and gfa ------
330C preconditioning: scale stress and strain so as to behave similarly to x,f -
331        do i = 1,3
332          do j = i,3
333            ndegi = ndegi + 1
334            gfa(ndegi) = -stress_dif(i,j)*volume/precon
335            gxa(ndegi) = strain(i,j) * precon
336            dxmaxn(ndegi) = ZmatMaxDisplLen
337          enddo
338        enddo
339
340C Check stress convergence
341        strtol = dabs(strtol)
342        do i = 1,3
343          do j = 1,3
344            relaxd = relaxd .and.
345     .        ( dabs(stress_dif(i,j)) .lt. abs(strtol) )
346          enddo
347        enddo
348
349      else   ! FIXED CELL
350
351C Set number of degrees of freedom & variables
352         relaxd = .true.
353        ndegi = 0
354          do i = 1,na
355            do k = 1,3
356              indi = 3*(i-1)+k
357              if (VaryZmat(indi)) then
358                ndegi = ndegi + 1
359                gxa(ndegi) = Zmat(indi)
360                vi = iZmattoVars(indi)
361                if (vi.eq.0) then
362                  force = ZmatForce(indi)
363                else
364                  force = ZmatForceVar(vi)
365                endif
366                gfa(ndegi) = force
367                if (ZmatType(indi).eq.1) then
368                  ftoln(ndegi) = ZmatForceTolAng
369                  dxmaxn(ndegi) = ZmatMaxDisplAng
370                else
371                  ftoln(ndegi) = ZmatForceTolLen
372                  dxmaxn(ndegi) = ZmatMaxDisplLen
373                endif
374                relaxd=relaxd.and.(dabs(force).lt.ftoln(ndegi))
375              endif
376            enddo
377          enddo
378
379      endif
380
381      if (relaxd) then
382
383         call local_dealloc()
384
385         return
386      end if
387
388      if (.not. broyden_is_setup(br)) then
389         call broyden_reset(br,ndeg,maxit,cycle_on_maxit,
390     $        jinv0,0.01_dp,dxmaxn)
391      endif
392
393      allocate(rnew(1:ndeg))
394
395      weight = 1.0_dp
396      call broyden_step(br,gxa,gfa,w=weight,newx=rnew)
397
398
399      ! Transform back
400      if ( varcel ) then
401
402        ! New cell
403        indi = ndeg-6
404        do i = 1,3
405          do j = i,3
406            indi = indi + 1
407            strain(i,j) = rnew(indi) / precon
408            strain(j,i) = strain(i,j)
409          enddo
410        enddo
411
412        cell = cellin + matmul(strain-1.0_dp,cellin)
413        if (constant_volume) then
414           new_volume = volcel(cell)
415           if (Node.eq.0) write(6,"(a,f12.4)")
416     $          "Volume before coercion: ",  new_volume/Ang**3
417           cell =  cell * (initial_volume/new_volume)**(1.0_dp/3.0_dp)
418        endif
419
420C Output fractional coordinates to cartesian Bohr, and copy to xa -----------
421        ndegi = 0
422        do ia=1,na
423            do k=1,3
424              indi = 3*(ia-1)+k
425              if (VaryZmat(indi)) then
426                ndegi = ndegi + 1
427                j = indi
428                do while (j.ne.0)
429                  if (ZmatType(j).gt.2) then
430                    Zmat(j) = 0.0_dp
431                    do n=1,3
432                      indi1 = 3*(ia-1)+n
433                      ! Assume all three coords of this atom
434                      ! are variables
435                      ndi = ndegi + n - k
436                      Zmat(j)=Zmat(j)+cell(k,n)*rnew(ndi)/modcel(n)
437                    enddo
438                  else
439                    Zmat(j) = rnew(ndegi)
440                  endif
441                  Zmat(j) = Zmat(j)*coeffA(j) + coeffB(j)
442                  j = iNextDept(j)
443                enddo
444              endif
445            enddo
446          enddo
447          call Zmat_to_Cartesian(xa)
448
449      else
450C Fixed cell
451C Return coordinates to correct arrays
452        ndegi = 0
453          do ia = 1,na
454            do k = 1,3
455              indi = 3*(ia-1)+k
456              if (VaryZmat(indi)) then
457                ndegi = ndegi + 1
458                j = indi
459                do while (j.ne.0)
460                  Zmat(j) = rnew(ndegi)*coeffA(j)+coeffB(j)
461                  j = iNextDept(j)
462                enddo
463              endif
464            enddo
465          enddo
466          call Zmat_to_Cartesian(xa)
467      endif
468
469!      print *, "zm_broyden: cell out : ", cell
470!      print *, "zm_broyden: strain out : ", strain
471
472      call local_dealloc()
473
474      contains
475
476      subroutine local_dealloc()
477
478C     Deallocate local memory
479      if ( allocated(dxmaxn) ) deallocate(dxmaxn)
480      if ( allocated(ftoln) ) deallocate(ftoln)
481      if ( allocated(gxa) ) deallocate(gxa)
482      if ( allocated(gfa) ) deallocate(gfa)
483      if ( allocated(rnew) ) deallocate(rnew)
484
485      end subroutine local_dealloc
486
487      end subroutine zm_broyden_optimizer
488
489      end module m_zm_broyden_optim
490