1c 2c $Id$ 3c 4 5* *************************** 6* * * 7* * md_energy * 8* * * 9* *************************** 10 real*8 function md_energy(E) 11 implicit none 12 real*8 E(*) 13 14#include "stdio.fh" 15#include "util.fh" 16#include "bafdecls.fh" 17#include "errquit.fh" 18 19 integer MASTER 20 parameter (MASTER=0) 21 22 logical value,oprint 23 integer taskid 24 integer i,j,ms 25 integer icount 26 real*8 EV,virial 27 real*8 cx,cy,cz 28 real*8 gx,gy,gz 29 real*8 en(2) 30 real*8 e_lj,e_q,e_spring,eion 31 integer rtdb 32 33* **** external functions **** 34 logical control_print 35 integer control_version 36 integer ion_nion,ion_katm,control_rtdb 37 real*8 ion_rion,ion_amass,ewald_e,ion_ion_e 38 external control_print 39 external control_version 40 external ion_nion,ion_katm,control_rtdb 41 external ion_rion,ion_amass,ewald_e,ion_ion_e 42 43 44* ***** QM/MM external functions **** 45 logical pspw_qmmm_found 46 real*8 pspw_qmmm_LJ_E,pspw_qmmm_mmq_Q_E,pspw_qmmm_spring_E 47 real*8 pspw_qmmm_LJ_Emix,pspw_qmmm_Q_Emix 48 external pspw_qmmm_found 49 external pspw_qmmm_LJ_E,pspw_qmmm_mmq_Q_E,pspw_qmmm_spring_E 50 external pspw_qmmm_LJ_Emix,pspw_qmmm_Q_Emix 51 52* ***** pspw_charge external functions **** 53 logical pspw_charge_found 54 real*8 pspw_charge_Energy_ion,pspw_charge_Energy_charge 55 external pspw_charge_found 56 external pspw_charge_Energy_ion,pspw_charge_Energy_charge 57 58 call Parallel_taskid(taskid) 59 oprint = ((taskid.eq.MASTER).and.control_print(print_medium)) 60 61* **** set the minimizer **** 62 !call dcopy(30,0.0d0,0,E,1) 63 64* **** generate phaze factors **** 65 !if (control_version().eq.3) call ewald_phafac() 66 67* ::::::::::: get energy ::::::::::::::::::::::: 68 eion = 0.0d0 69 if (control_version().eq.3) eion = ewald_e() 70 if (control_version().eq.4) eion = ion_ion_e() 71 72 73 E(1) = 0.0d0 74 E(2) = eion 75 E(3) = 0.0d0 76 E(4) = 0.0d0 77 E(5) = eion 78 E(6) = 0.0d0 79 E(7) = 0.0d0 80 E(8) = 0.0d0 81 E(9) = 0.0d0 82 E(10) = 0.0d0 83 if (pspw_qmmm_found()) then 84 e_lj = pspw_qmmm_LJ_E() 85 e_q = pspw_qmmm_mmq_Q_E() 86 e_spring = pspw_qmmm_spring_E() 87 E(2) = E(2) + e_lj + e_q + e_spring 88 E(11) = e_lj 89 E(12) = e_q 90 E(13) = e_spring 91 E(14) = pspw_qmmm_LJ_Emix() 92 E(14) = E(14) + pspw_qmmm_Q_Emix() 93 !E(23) = E(23) + E(14) 94 !E(24) = E(24) + E(14)*E(14) 95 96 end if 97 98 !E(25) = E(25) + E(1) 99 !E(26) = E(26) + E(1)*E(1) 100 101 102 103* **** get pspw_charge energies **** 104c if (pspw_charge_found()) then 105c E(19) = 0.0d0 106c E(20) = pspw_charge_Energy_ion() 107c E(21) = pspw_charge_Energy_charge() 108c E(1) = E(1) + E(20) + E(21) 109c end if 110 111 112c*::::::::::::::::: report summary of results ::::::::::::::::::::::: 113c if (oprint) then 114c write(luout,1304) 115c write(luout,1410) 116c 117c write(luout,*) 118c write(luout,1430) E(1),E(1)/ion_nion() 119c 120c 121c if (pspw_charge_found()) then 122c write(luout,1431) 123c write(luout,1432) 124c write(luout,1433) (E(1)-E(19)-E(20)-E(21)), 125c > (E(1)-E(19)-E(20)-E(21))/ion_nion() 126c end if 127c 128c write(luout,1470) E(5),E(5)/ion_nion() 129c if (pspw_qmmm_found()) then 130c write(luout,1700) 131c write(luout,1701) 132c write(luout,1702) E(11) 133c write(luout,1703) E(12) 134c write(luout,1704) E(13) 135c end if 136c if (pspw_charge_found()) then 137c write(luout,1800) 138c write(luout,1801) 139c write(luout,1805) E(19)+E(20)+E(21) 140c write(luout,1802) E(19) 141c write(luout,1803) E(20) 142c write(luout,1804) E(21) 143c end if 144c 145c end if 146c call ecce_print1 ('total energy', mt_dbl, E(1), 1) 147c call ecce_print1 ('nuclear repulsion energy', mt_dbl, E(5), 1) 148 149 md_energy = E(2) 150 return 151 152 1190 FORMAT(5X, I4, A5 ,' (',3F11.5,' ) - atomic mass= ',F6.3,' ') 153 1200 FORMAT(5X,' G.C. ',' (',3F11.5,' )') 154 1210 FORMAT(5X,' C.O.M.',' (',3F11.5,' )') 155 1300 FORMAT(//'======================') 156 1301 FORMAT(//'== Energy Calculation ==') 157 1302 FORMAT( '======================') 158 1304 FORMAT(/) 159 1310 FORMAT(I8,E20.10,3E15.5) 160 1320 FORMAT(' number of electrons: spin up=',F11.5,' down=',F11.5,A) 161 1350 FORMAT(/' orthonormality') 162 1360 FORMAT(I3,2I3,E18.7) 163 1370 FORMAT(I3) 164 1380 FORMAT(' ''',a,'''',I4) 165 1390 FORMAT(I3) 166 1400 FORMAT(I3,3E18.8/3X,3E18.8) 167c1410 FORMAT(10X,'============= summary of results =================') 168 1410 FORMAT('== Summary Of Results ==') 169 1420 FORMAT( ' final position of ions:') 170 1430 FORMAT(/' total energy :',E19.10,' (',E15.5,'/ion)') 171 1431 FORMAT(/' QM Energies') 172 1432 FORMAT( '------------') 173 1433 FORMAT( ' total QM energy :',E19.10,' (',E15.5,'/ion)') 174 1440 FORMAT( ' total orbital energy:',E19.10,' (',E15.5,'/electron)') 175 1450 FORMAT( ' hartree energy :',E19.10,' (',E15.5,'/electron)') 176 1455 FORMAT( ' SIC-hartree energy :',E19.10,' (',E15.5,'/electron)') 177 1456 FORMAT( ' SIC-exc-corr energy :',E19.10,' (',E15.5,'/electron)') 178 1457 FORMAT( ' HF exchange energy :',E19.10,' (',E15.5,'/electron)') 179 1460 FORMAT( ' exc-corr energy :',E19.10,' (',E15.5,'/electron)') 180 1470 FORMAT( ' ion-ion energy :',E19.10,' (',E15.5,'/ion)') 181 1471 FORMAT( ' smearing energy :',E19.10,' (',E15.5,'/electron)') 182 1480 FORMAT(/' kinetic (planewave) :',E19.10,' (',E15.5,'/electron)') 183 1490 FORMAT( ' V_local (planewave) :',E19.10,' (',E15.5,'/electron)') 184 1491 FORMAT( ' Vl+Vqm/mm :',E19.10,' (',E15.5,'/electron)') 185 1495 FORMAT( ' V_nl (planewave) :',E19.10,' (',E15.5,'/electron)') 186 1496 FORMAT( ' V_Coul (planewave) :',E19.10,' (',E15.5,'/electron)') 187 1497 FORMAT( ' V_xc. (planewave) :',E19.10,' (',E15.5,'/electron)') 188 1498 FORMAT( ' Virial Coefficient :',E19.10) 189 1499 FORMAT( ' K.S. SIC-hartree energy :',E19.10, 190 > ' (',E15.5,'/electron)') 191 1501 FORMAT( ' K.S. SIC-exc-corr energy :',E19.10, 192 > ' (',E15.5,'/electron)') 193 1502 FORMAT( ' K.S. HFX energy :',E19.10, 194 > ' (',E15.5,'/electron)') 195 1500 FORMAT(/' orbital energies:') 196 1507 FORMAT(/' Fermi energy =',2(E18.7,' (',F8.3,'eV)')) 197 1510 FORMAT(2(E18.7,' (',F8.3,'eV)')) 198 1511 FORMAT(2(E18.7,' (',F8.3,'eV) occ=',F5.3)) 199 1512 FORMAT(2(E18.7,' (',F8.3,'eV)',A4)) 200 1513 FORMAT(2(E18.7,' (',F8.3,'eV)',A4,' occ=',F5.3)) 201 202 203 1700 FORMAT(/' LJ/residualCoulomb/vib/CAV Energies') 204 1701 FORMAT( ' -----------------------------------') 205 1702 FORMAT( ' LJ energy :',E19.10) 206 1703 FORMAT( ' Residual Coulomb energy :',E19.10) 207 1704 FORMAT( ' MM Vibrational energy :',E19.10) 208 1705 FORMAT( ' MM Vibration energy :',E19.10) 209 1706 FORMAT( ' (QM+MM)/Cavity energy :',E19.10) 210 1707 FORMAT( ' - MM Charge Field/QM Electron :',E19.10) 211 1708 FORMAT( ' - MM Charge Field/QM Ion :',E19.10) 212 1709 FORMAT( ' - MM LJ/QM LJ :',E19.10) 213 1710 FORMAT( ' - MM Charge Field/MM Charge Field:',E19.10) 214 1711 FORMAT( ' - MM LJ/MM LJ :',E19.10) 215 216 1800 FORMAT(/' Charge Field Energies') 217 1801 FORMAT( ' ---------------------') 218 1802 FORMAT( ' - Charge Field/Electron :',E19.10) 219 1803 FORMAT( ' - Charge Field/Ion :',E19.10) 220 1804 FORMAT( ' - Charge Field/Charge Field:',E19.10) 221 1805 FORMAT( ' Charge Field Energy :',E19.10) 222 end 223 224 225 226* ******************************* 227* * * 228* * md_energy_gradient * 229* * * 230* ******************************* 231 subroutine md_energy_gradient(G1) 232 implicit none 233 real*8 G1(3,*) 234 235#include "stdio.fh" 236#include "util.fh" 237 238 logical allow_translation,lprint,mprint 239 integer MASTER 240 parameter (MASTER=0) 241 integer i,k,taskid,nion,nion1 242 integer i1 243 real*8 GG,fmax,fatom 244 real*8 fmx,fmy,fmz 245 real*8 fmx2,fmy2,fmz2 246 247* **** external functions **** 248 logical psp_semicore,pspw_charge_found,pspw_qmmm_found 249 logical control_allow_translation,ion_q_FixIon,control_print 250 character*4 ion_aname,pspw_charge_aname 251 integer ion_katm,ion_nion,control_version 252 integer pspw_charge_nion 253 real*8 ion_rion,pspw_charge_rion 254 real*8 pspw_charge_charge 255 external psp_semicore,pspw_charge_found,pspw_qmmm_found 256 external control_allow_translation,ion_q_FixIon,control_print 257 external ion_aname,pspw_charge_aname 258 external ion_katm,ion_nion,control_version 259 external pspw_charge_nion 260 external ion_rion,pspw_charge_rion 261 external pspw_charge_charge 262 logical pspw_bqext 263 external pspw_bqext 264 !integer ion_nion_mm,ion_nion_qm 265 !external ion_nion_mm,ion_nion_qm 266 267 allow_translation = control_allow_translation() 268 nion = ion_nion() 269 if (pspw_charge_found().and. 270 > (.not.pspw_bqext())) nion = nion + pspw_charge_nion() 271 272* **** generate phaze factors **** 273 if (control_version().eq.3) call ewald_phafac() 274 275 call dcopy(3*nion,0.0d0,0,G1,1) 276 if (control_version().eq.3) call ewald_f(G1) 277 if (control_version().eq.4) call ion_ion_f(G1) 278 if (pspw_qmmm_found()) call pspw_qmmm_mmq_fion(G1) 279 280c if (pspw_charge_found()) then 281c if(pspw_bqext()) then 282c call pspw_charge_charge_Fion(G1) 283c else 284c nion1 = ion_nion() 285c call pspw_charge_Fion_Fcharge(G1,G1(1,nion1+1)) 286c call pspw_charge_Fcharge(G1(1,nion1+1)) 287c end if 288c end if 289 290* **** remove ion forces using ion_FixIon **** 291 call ion_FixIon(G1) 292 293 if (.not.allow_translation) then 294 call center_F_mass(G1,fmx,fmy,fmz) 295 do i=1,nion 296 G1(1,i) = G1(1,i) - fmx 297 G1(2,i) = G1(2,i) - fmy 298 G1(3,i) = G1(3,i) - fmz 299 end do 300 end if 301 call center_F_mass(G1,fmx2,fmy2,fmz2) 302 303 GG = 0.0d0 304 fmax = 0.0d0 305 do i=1,nion 306 GG = GG + G1(1,i)**2 + G1(2,i)**2 + G1(3,i)**2 307 fatom = dsqrt(G1(1,i)**2+G1(2,i)**2 +G1(3,i)**2) 308 if (fatom.gt.fmax) fmax = fatom 309 end do 310 311 call Parallel_taskid(taskid) 312 mprint = ((taskid.eq.MASTER).and.control_print(print_high)) 313 lprint = ((taskid.eq.MASTER).and.control_print(print_medium)) 314 315 if (taskid.eq.MASTER) then 316 if (mprint) write(luout,1301) 317 if (lprint) then 318 write(luout,1304) 319 if (.not.allow_translation) write(luout,1400) fmx,fmy,fmz 320 write(luout,1304) 321 write(luout,1410) 322 end if 323 324* **** print out positions *** 325 if (mprint) then 326 write(luout,1420) 327 do I=1,ion_nion() 328 if (ion_q_FixIon(I)) then 329 write(6,1191) I,ion_aname(I),(ion_rion(K,I),K=1,3) 330 else 331 write(6,1190) I,ion_aname(I),(ion_rion(K,I),K=1,3) 332 end if 333 end do 334 335c* **** print out charge positions *** 336c if (pspw_charge_found().and.(.not.pspw_bqext())) then 337c do i=1,pspw_charge_nion() 338c i1 = ion_nion() + i 339c if (ion_q_FixIon(i1)) then 340c write(luout,1193) i1,pspw_charge_aname(i), 341c > (pspw_charge_rion(K,i),K=1,3), 342c > pspw_charge_charge(i) 343c else 344c write(luout,1192) i1,pspw_charge_aname(i), 345c > (pspw_charge_rion(K,i),K=1,3), 346c > pspw_charge_charge(i) 347c end if 348c end do 349c end if 350 351 end if 352 353c* **** print out forces *** 354 if (lprint) then 355 write(luout,1421) 356 write(luout,1190)(i,ion_aname(I), 357 > (G1(K,I),K=1,3),I=1,ion_nion()) 358 359c* **** print out charge forces *** 360c if (pspw_charge_found().and.(.not.pspw_bqext())) then 361c do i=1,pspw_charge_nion() 362c i1 = ion_nion() + i 363c write(luout,1190) i1,pspw_charge_aname(i), 364c > (G1(K,i1),K=1,3) 365c end do 366c end if 367c 368c write(luout,1210) fmx2,fmy2,fmz2 369c write(luout,1425) 370c write(luout,1426) dsqrt(GG), 371c > dsqrt(GG)/dble(nion), 372c > fmax,fmax*(27.2116d0/0.529177d0) 373 374 end if 375 end if 376 377c call dscal(3*nion,(-1.0d0),G1,1) 378 379 return 380 1190 FORMAT(5X, I4, A5, ' (',3F11.5,' )') 381 1191 FORMAT(5X, I4, A5, ' (',3F11.5,' ) - fixed') 382 1192 FORMAT(5X, I4, A5, ' (',3F11.5,' ) q=',F8.3) 383 1193 FORMAT(5X, I4, A5, ' (',3F11.5,' ) q=',F8.3,' - fixed') 384 1210 FORMAT(5X,' C.O.M.',' (',3F11.5,' )') 385 1300 FORMAT(//'========================') 386 1301 FORMAT(//'== Gradient Calculation ==') 387 1302 FORMAT( '========================') 388 1304 FORMAT(/) 389 1400 FORMAT('Translation force removed: (',3F11.5,')') 390 1410 FORMAT(10X,'============= Ion Gradients =================') 391 1425 FORMAT(10X,'===============================================') 392 1426 FORMAT(10X,'|F| =',E15.6, 393 > /10x,'|F|/nion =',E15.6, 394 > /10x,'max|Fatom|=',E15.6,1x,'(',F8.3,'eV/Angstrom)'//) 395 1420 FORMAT( ' Ion Positions:') 396 1421 FORMAT( ' Ion Forces:') 397 end 398 399