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