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