1c
2c     $Id$
3c
4
5*  **************************************************************
6*  *                MPI cgminimize2 routine                     *
7*  *               (Fletcher-Reeves' steps)                     *
8*  *                                                            *
9*  *  This is a developing Stiefel conjugate gradient minimizer *
10*  *  written for NWChem.                                       *
11*  *                                                            *
12*  *                                                            *
13*  **************************************************************
14
15      subroutine cgminimize2(E,deltae,deltac,current_iteration)
16      implicit none
17      real*8     E(*)
18      real*8     deltae,deltac
19      integer    current_iteration
20
21#include "bafdecls.fh"
22#include "errquit.fh"
23
24*     **** local variables ****
25
26      real*8  deltat_min
27      parameter (deltat_min=1.0d-3)
28
29      integer H0(2),G1(2)
30      real*8  E0,dE0
31
32      real*8     sum0(2),sum1(2),scale,tole,tolc
33      real*8     ehartree,eorbit,exc,pxc,eion
34      real*8     Enew,Eold,Estart
35      common / cgsd_block / Enew,Eold,Estart
36
37      integer it,it_in,ms,ispin,ne(2)
38      real*8 tmin,deltat
39      real*8 max_sigma
40
41      logical value
42      integer neall,npack1
43      !real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
44      !real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
45      real*8 e_lj,e_q,e_spring
46      real*8 ehsic,exsic,phsic,pxsic,ehfx,phfx
47
48
49*     **** external functions ****
50      integer  control_it_in,psi_neq,control_version,control_ispin
51      real*8   control_tole,control_tolc
52      real*8   psi_geodesic2_energy
53      real*8   psi_geodesic2_denergy
54      real*8   rho_error
55      real*8   dng_1ehartree
56      real*8   psi_1ke
57      real*8   psi_1vl,psi_1v_field,dng_1vl_mm
58      real*8   psi_1vnl
59      real*8   rho_1exc
60      real*8   rho_1pxc
61      real*8   ewald_e,ion_ion_e
62      real*8   psi_1eorbit
63      real*8   linesearch
64
65      external control_it_in,psi_neq,control_version,control_ispin
66      external control_tole,control_tolc
67      external psi_geodesic2_energy
68      external psi_geodesic2_denergy
69      external rho_error
70      external dng_1ehartree
71      external psi_1ke
72      external psi_1vl,psi_1v_field,dng_1vl_mm
73      external psi_1vnl
74      external rho_1exc
75      external rho_1pxc
76      external ewald_e,ion_ion_e
77      external psi_1eorbit
78      external linesearch
79
80*     ***** QM/MM external functions ****
81      logical  pspw_qmmm_found
82      real*8   pspw_qmmm_LJ_E,pspw_qmmm_Q_E,pspw_qmmm_spring_E
83      real*8   pspw_qmmm_LJ_Emix,pspw_qmmm_Q_Emix
84      external pspw_qmmm_found
85      external pspw_qmmm_LJ_E,pspw_qmmm_Q_E,pspw_qmmm_spring_E
86      external pspw_qmmm_LJ_Emix,pspw_qmmm_Q_Emix
87
88*     ***** pspw_charge external functions ****
89      logical  pspw_charge_found
90      real*8   pspw_charge_Energy_ion,pspw_charge_Energy_charge
91      external pspw_charge_found
92      external pspw_charge_Energy_ion,pspw_charge_Energy_charge
93
94      logical  control_precondition,pspw_SIC,pspw_HFX,psp_U_psputerm
95      external control_precondition,pspw_SIC,pspw_HFX,psp_U_psputerm
96      logical  meta_found,ion_disp_on
97      external meta_found,ion_disp_on
98      real*8   ion_disp_energy,psi_smearcorrection
99      external ion_disp_energy,psi_smearcorrection
100      logical  control_fractional
101      external control_fractional
102
103
104      call Pack_npack(1,npack1)
105      ispin = control_ispin()
106      ne(1) = psi_neq(1)
107      ne(2) = psi_neq(2)
108      neall = ne(1)+ne(2)
109
110*     **** check and fix orthogonality ****
111      call psi_1ortho_check_fix()
112
113
114*     **** allocate H0 and G1 ****
115      value = BA_alloc_get(mt_dcpl,npack1*neall,
116     >                     'H0',H0(2),H0(1))
117      value = value.and.
118     >        BA_alloc_get(mt_dcpl,npack1*neall,
119     >                     'G1',G1(2),G1(1))
120      if (.not. value) call errquit('cgminimize2:out of heap memory',0,
121     &       MA_ERR)
122      call dcopy(2*npack1*neall,0.0d0,0,dcpl_mb(G1(1)),1)
123
124      Estart = Enew
125
126*     ***** get the initial gradient and direction ****
127      call psi_1get_TSgradient(dcpl_mb(G1(1)),E0)
128
129
130      do ms=1,ispin
131        call Grsm_gg_trace(npack1,ne(ms),
132     >                     dcpl_mb(G1(1)+(ms-1)*npack1*ne(1)),
133     >                     dcpl_mb(G1(1)+(ms-1)*npack1*ne(1)),
134     >                     sum1(ms))
135        call D1dB_SumAll(sum1(ms))
136      end do
137
138      call Grsm_gg_Copy(npack1,neall,
139     >                   dcpl_mb(G1(1)),
140     >                   dcpl_mb(H0(1)))
141
142*     ******************************************
143*     ****                                  ****
144*     **** Start of conjugate gradient loop ****
145*     ****                                  ****
146*     ******************************************
147      it_in = control_it_in()
148      tole  = control_tole()
149      tolc  = control_tolc()
150      tmin  = deltat_min
151      do it=2,it_in
152
153*        **** initialize the geoedesic line data structure ****
154         call psi_1geodesic2_start(dcpl_mb(H0(1)),max_sigma,dE0)
155
156
157*        ******* line search *********
158         if (tmin.gt.deltat_min) then
159            deltat = tmin
160         else
161            deltat = deltat_min
162         end if
163         Enew = linesearch(0.0d0,E0,dE0,deltat,
164     >                        psi_geodesic2_energy,
165     >                        psi_geodesic2_denergy,
166     >                        0.5d0,tmin,deltae,2)
167         call psi_geodesic2_final(tmin)
168         deltac = rho_error()
169
170*        **** exit loop early ****
171         if ((dabs(deltae).lt.tole).and.(deltac.lt.tolc)) then
172            go to 30
173         end if
174
175
176*        **** transport the previous search directions ****
177         call psi_1geodesic2_transport(tmin,dcpl_mb(H0(1)))
178
179*        **** make psi1 <--- psi2(tmin) ****
180         call psi_2to1()
181
182*        **** get the new gradient - also updates densities****
183         call psi_1get_TSgradient(dcpl_mb(G1(1)),E0)
184
185
186         do ms=1,ispin
187           sum0(ms)  = sum1(ms)
188           call Grsm_gg_trace(npack1,ne(ms),
189     >                        dcpl_mb(G1(1)+(ms-1)*npack1*ne(1)),
190     >                        dcpl_mb(G1(1)+(ms-1)*npack1*ne(1)),
191     >                        sum1(ms))
192           call D1dB_SumAll(sum1(ms))
193         end do
194
195*        **** the new direction using Fletcher-Reeves ****
196         if ( (dabs(deltae).le.(1.0d-2)).and.
197     >        (tmin.gt.deltat_min)) then
198
199           do ms=1,ispin
200             if (sum0(ms).gt.1.0d-6) then
201                scale = sum1(ms)/sum0(ms)  ! Fletcher-Reeves
202             else
203                scale = 0.0d0
204             end if
205
206c             call Grsm_gg_dScale(npack1,ne(ms),scale,
207c     >                         dcpl_mb(H0(1)+(ms-1)*npack1*ne(1)),
208c     >                         dcpl_mb(H0(1)+(ms-1)*npack1*ne(1)))
209             call Grsm_gg_dScale1(npack1,ne(ms),scale,
210     >                         dcpl_mb(H0(1)+(ms-1)*npack1*ne(1)))
211           end do
212c           call Grsm_ggg_Sum(npack1,neall,
213c     >                       dcpl_mb(G1(1)),
214c     >                       dcpl_mb(H0(1)),
215c     >                       dcpl_mb(H0(1)))
216           call Grsm_ggg_Sum2(npack1,neall,
217     >                       dcpl_mb(G1(1)),
218     >                       dcpl_mb(H0(1)))
219
220
221*        **** the new direction using steepest-descent ****
222         else
223              call Grsm_gg_Copy(npack1,neall,
224     >                          dcpl_mb(G1(1)),
225     >                          dcpl_mb(H0(1)))
226         end if
227
228      end do
229
230
231*     **** initialize the geoedesic line data structure ****
232      call psi_1geodesic2_start(dcpl_mb(H0(1)),max_sigma,dE0)
233
234
235*     ******* line search *********
236      if (tmin.gt.deltat_min) then
237         deltat = tmin
238      else
239         deltat = deltat_min
240      end if
241
242      Enew = linesearch(0.0d0,E0,dE0,deltat,
243     >                        psi_geodesic2_energy,
244     >                        psi_geodesic2_denergy,
245     >                        0.5d0,tmin,deltae,2)
246
247      call psi_geodesic2_final(tmin)
248      deltac = rho_error()
249
250 30   call psi_2to1()
251      call psi_1toelectron()
252c      call psi_check()
253
254      eion = 0.0d0
255      if (control_version().eq.3) eion = ewald_e()
256      if (control_version().eq.4) eion = ion_ion_e()
257
258      eorbit   = psi_1eorbit()
259      ehartree = dng_1ehartree()
260      exc      = rho_1exc()
261      pxc      = rho_1pxc()
262
263      E(1)  = Enew + eion
264      E(2)  = eorbit
265      E(3)  = ehartree
266      E(4)  = exc
267      E(5)  = eion
268      E(6)  = psi_1ke()
269      E(7)  = psi_1vl()
270      E(8)  = psi_1vnl()
271      E(9)  = 2.0d0*ehartree
272      E(10) = pxc
273      if (pspw_qmmm_found()) then
274         e_lj     = pspw_qmmm_LJ_E()
275         e_q      = pspw_qmmm_Q_E()
276         e_spring = pspw_qmmm_spring_E()
277         E(1)  = E(1) + e_lj + e_q + e_spring
278
279         E(11) = e_lj
280         E(12) = e_q
281         E(13) = e_spring
282
283*        **** Eqm-mm energy ***
284         E(14) = pspw_qmmm_LJ_Emix()
285         E(14) = E(14) + pspw_qmmm_Q_Emix()
286         E(14) = E(14) + dng_1vl_mm()
287      end if
288
289
290
291*     **** get pspw_charge  energies ****
292      if (pspw_charge_found()) then
293         E(19)  = psi_1v_field()
294         E(20)  = pspw_charge_Energy_ion()
295         E(21)  = pspw_charge_Energy_charge()
296         E(1)   = E(1) + E(20) + E(21)
297      end if
298
299
300*     **** SIC corrections ****
301      if (pspw_SIC()) then
302         call electron_SIC_energies(ehsic,phsic,exsic,pxsic)
303         E(22) = ehsic
304         E(23) = exsic
305         E(24) = phsic
306         E(25) = pxsic
307      end if
308
309*     **** HFX terms ****
310      if (pspw_HFX()) then
311         call electron_HFX_energies(ehfx,phfx)
312         E(26) = ehfx
313         E(27) = phfx
314      end if
315
316*     **** DFT+U terms ****
317      if (psp_U_psputerm()) then
318         call electron_U_energies(ehfx,phfx)
319         E(29) =  ehfx
320         E(30) =  phfx
321      end if
322
323*     **** Metadynamics potential terms ****
324      if (meta_found()) then
325         call electron_meta_energies(ehfx,phfx)
326         E(31) =  ehfx
327         E(32) =  phfx
328      end if
329
330      if (control_fractional()) then
331        E(28) = psi_smearcorrection()
332        E(1)  = E(1) + E(28)
333      end if
334
335*     **** Dispersion energy ****
336      if (ion_disp_on()) then
337         E(33) = ion_disp_energy()
338         E(1)  = E(1) + E(33)
339      end if
340
341      value = BA_free_heap(G1(2))
342      value = value.and.BA_free_heap(H0(2))
343      if (.not. value)
344     >  call errquit('cgminimize2:error freeing heap memory',2, MA_ERR)
345
346
347      return
348      end
349
350
351