1*
2* $Id$
3*
4***********************************************************************
5*                     band_minimizer								  *
6*                                                                     *
7*     This is a developing band structure parallel code for NWCHEM    *
8*       + tcgmsg message passing library used                         *
9*       + my own slap-decomposed parallel 3d-FFT(real->complex) used  *
10*                                                                     *
11*                                                                     *
12***********************************************************************
13
14      logical function band_minimizer(rtdb,flag)
15      implicit none
16      integer rtdb
17      integer flag
18
19#include "global.fh"
20#include "bafdecls.fh"
21#include "inp.fh"
22#include "btdb.fh"
23#include "stdio.fh"
24#include "util.fh"
25#include "errquit.fh"
26
27*     **** parallel variables ****
28      integer  taskid,np,np_i,np_j,np_k
29      integer  MASTER
30      parameter(MASTER=0)
31
32*     **** timing variables ****
33      real*8   cpu1,cpu2,cpu3,cpu4
34      real*8   t1,t2,t3,t4,av
35
36*     **** lattice variables ****
37      integer ngrid(3),nwave,nfft3d
38      real*8  a,b,c,alpha,beta,gamma
39
40*     ***** energy variables ****
41      logical spin_orbit
42      integer ispin,ne(2)
43      real*8  E(10)
44      real*8  dipole(3)
45      real*8  stress(3,3),lstress(6)
46
47*     **** gradient variables ****
48      integer fion(2)
49
50*     **** error variables ****
51      logical value
52      integer ierr
53
54*     **** local variables ****
55      logical newpsi,lprint,mprint,hprint
56      real*8  gx,gy,gz,cx,cy,cz
57      real*8  EV,pi,weight
58      real*8  icharge,en(2)
59      real*8  f0,f1,f2,f3,f4,f5,f6
60      integer if1,if2
61      integer i,k,ia,nion,vers,nbrillioun,nb,ms
62      integer mapping,minimizer
63
64*     **** external functions ****
65      real*8      lattice_omega,lattice_unita,lattice_ecut,lattice_wcut
66      real*8      lattice_unitg,ion_amass,ion_TotalCharge
67      character*4 ion_aname
68      external    lattice_omega,lattice_unita,lattice_ecut,lattice_wcut
69      external    lattice_unitg,ion_amass,ion_TotalCharge
70      external    ion_aname
71
72
73      real*8   control_tole,control_tolc,control_tolr,ion_rion
74      external control_tole,control_tolc,control_tolr,ion_rion
75      real*8   control_time_step,control_fake_mass
76      external control_time_step,control_fake_mass
77      logical  control_read,control_move,ion_init,ion_q_FixIon
78      external control_read,control_move,ion_init,ion_q_FixIon
79      logical control_spin_orbit
80      external control_spin_orbit
81      integer  control_it_in,control_it_out,control_gga,control_version
82      integer  control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm
83      integer  ion_nkatm
84      external control_it_in,control_it_out,control_gga,control_version
85      external control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm
86      external ion_nkatm
87
88      character*12 control_boundry
89      external     control_boundry
90
91      logical  brillioun_print,control_Mulliken
92      integer  brillioun_nbrillioun,brillioun_nbrillq
93      integer  brillioun_nbrillioun0
94      real*8   brillioun_weight_brdcst,brillioun_ks_brdcst
95      real*8   brillioun_k_brdcst,brillioun_weight
96      external brillioun_print,control_Mulliken
97      external brillioun_nbrillioun,brillioun_nbrillq
98      external brillioun_nbrillioun0
99      external brillioun_weight_brdcst,brillioun_ks_brdcst
100      external brillioun_k_brdcst,brillioun_weight
101      integer  c_electron_count,linesearch_count
102      external c_electron_count,linesearch_count
103
104      real*8   nwpw_timing
105      external nwpw_timing
106      integer  Cram_nwave_all_brdcst,Cram_nwave_brdcst
107      external Cram_nwave_all_brdcst,Cram_nwave_brdcst
108
109      integer  ewald_ncut
110      real*8   ewald_rcut,ewald_mandelung,ewald_e
111      external ewald_ncut
112      external ewald_rcut,ewald_mandelung,ewald_e
113      logical  cpsp_semicore,psi_filefind,cpsi_initialize,cpsi_finalize
114      external cpsp_semicore,psi_filefind,cpsi_initialize,cpsi_finalize
115      real*8   c_cgsd_energy,c_cgsd_noit_energy
116      external c_cgsd_energy,c_cgsd_noit_energy
117      integer  cpsp_lmax,cpsp_locp,cpsp_lmmax
118      external cpsp_lmax,cpsp_locp,cpsp_lmmax
119      integer  cpsp_nprj,cpsp_psp_type
120      external cpsp_nprj,cpsp_psp_type
121      real*8   cpsp_rcore,cpsp_rc,cpsp_ncore,cpsp_zv
122      external cpsp_rcore,cpsp_rc,cpsp_ncore,cpsp_zv
123      character*4 ion_atom
124      external    ion_atom
125
126      character*255 cpsp_comment,comment
127      external      cpsp_comment
128
129      integer  cpsi_ispin,cpsi_ne,psi_get_version
130      external cpsi_ispin,cpsi_ne,psi_get_version
131      logical  pspw_reformat_c_wvfnc,control_check_charge_multiplicity
132      external pspw_reformat_c_wvfnc,control_check_charge_multiplicity
133      integer  control_mapping,control_minimizer,control_np_dimensions
134      external control_mapping,control_minimizer,control_np_dimensions
135      integer  control_scf_algorithm
136      external control_scf_algorithm
137      real*8   control_ks_alpha,control_kerker_g0
138      external control_ks_alpha,control_kerker_g0
139      integer  control_ks_maxit_orb,control_ks_maxit_orbs
140      external control_ks_maxit_orb,control_ks_maxit_orbs
141      real*8   cpsi_occupation,control_fractional_temperature
142      external cpsi_occupation,control_fractional_temperature
143      logical  control_fractional,control_print
144      external control_fractional,control_print
145      integer  control_fractional_smeartype
146      external control_fractional_smeartype
147      real*8   control_fractional_kT,control_fractional_alpha
148      external control_fractional_kT,control_fractional_alpha
149
150c      character*255 cpsp_comment,comment
151c      external      cpsp_comment
152
153      integer  ion_nconstraints,ion_ndof
154      external ion_nconstraints,ion_ndof
155      logical  control_smooth_cutoff
156      external control_smooth_cutoff
157      real*8   control_smooth_cutoff_values
158      external control_smooth_cutoff_values
159
160
161*****************************|  PROLOGUE  |****************************
162
163      value = .true.
164      pi = 4.0d0*datan(1.0d0)
165
166      call nwpw_timing_init()
167      call dcopy(10,0.0d0,0,E,1)
168
169
170*     **** get parallel variables ****
171      call Parallel_Init()
172      call Parallel_np(np)
173      call Parallel_taskid(taskid)
174
175      value = control_read(5,rtdb)
176      if (.not. value)
177     > call errquit('error reading control',0, DISK_ERR)
178
179      lprint = ((taskid.eq.MASTER).and.(control_print(print_low)))
180      mprint = ((taskid.eq.MASTER).and.(control_print(print_medium)))
181      hprint = ((taskid.eq.MASTER).and.(control_print(print_high)))
182
183      if (taskid.eq.MASTER) call current_second(cpu1)
184*     ***** print out header ****
185      if (mprint) then
186         write(luout,1000)
187         write(luout,1010)
188         write(luout,1020)
189         write(luout,1010)
190         write(luout,1030)
191         write(luout,1010)
192         write(luout,1035)
193         write(luout,1010)
194         write(luout,1040)
195         write(luout,1010)
196         write(luout,1041)
197         write(luout,1042)
198         write(luout,1043)
199         write(luout,1044)
200         write(luout,1010)
201         write(luout,1000)
202         call nwpw_message(1)
203         write(luout,1110)
204         call flush(luout)
205      end if
206
207
208      call Parallel3d_Init(control_np_dimensions(2),
209     >                     control_np_dimensions(3))
210      call Parallel3d_np_i(np_i)
211      call Parallel3d_np_j(np_j)
212      call Parallel3d_np_k(np_k)
213
214      ngrid(1) = control_ngrid(1)
215      ngrid(2) = control_ngrid(2)
216      ngrid(3) = control_ngrid(3)
217      nwave = 0
218      minimizer = control_minimizer()
219      mapping   = control_mapping()
220      ierr = 0
221
222*     **** initialize C3dB data structure ****
223      call C3dB_Init(1,ngrid(1),ngrid(2),ngrid(3),mapping)
224      call C3dB_nfft3d(1,nfft3d)
225
226      call cpsi_data_init(20)
227
228*     **** read ions ****
229      value = ion_init(rtdb)
230      call center_geom(cx,cy,cz)
231      call center_mass(gx,gy,gz)
232
233*     **** initialize lattice data structure ****
234      call lattice_init()
235      call c_G_init()
236
237*     **** nitialize brillion ***
238      call brillioun_init()
239      call Cram_Init()
240      call C3dB_pfft_init()
241
242*     **** initialize D3dB data structure and mask for GGA ****
243      if ((control_gga().ge.10).and.(control_gga().le.200)) THEN
244      call D3dB_Init(1,ngrid(1),ngrid(2),ngrid(3),mapping)
245      call G_init()
246      call mask_init()
247      end if
248
249c*     **** read ions ****
250c      value = ion_init(rtdb)
251c      call center_geom(cx,cy,cz)
252c      call center_mass(gx,gy,gz)
253
254*     **** allocate psp data structure and read in psedupotentials into it ****
255      call cpsp_init()
256      call cpsp_readall()
257      if (cpsp_semicore(0)) call c_semicore_check()
258
259
260*     **** initialize ke,and coulomb data structures ****
261      call cstrfac_init()
262      call cke_init()
263      call c_coulomb_init()
264      call ewald_init()
265
266*     **** set up phase factors at the current geometry  ****
267      call cphafac()
268      call cphafac_k()
269      call ewald_phafac()
270
271*     **** read in wavefunctions and initialize psi ****
272      if (.not.control_check_charge_multiplicity()) then
273         call cpsi_new()
274         newpsi = .true.
275      else
276         newpsi = .false.
277*        **** convert from pspw format to band format ****
278         vers = psi_get_version()
279         if ((vers.eq.3).or.(vers.eq.4)) then
280            newpsi = .true.
281            value = btdb_parallel(.false.)
282            if (taskid.eq.MASTER) then
283              value= pspw_reformat_c_wvfnc(1)
284            end if
285            value = btdb_parallel(.true.)
286         end if
287      end if
288
289      call psi_get_ne(ispin,ne)
290      if (ispin.eq.3) then
291         spin_orbit = .true.
292         ispin=2
293      else
294         spin_orbit = .false.
295      end if
296      nbrillioun = brillioun_nbrillioun()
297      call Pneb_init(ispin,ne,nbrillioun,spin_orbit)
298      value = cpsi_initialize(.true.)
299      newpsi = newpsi.or.value
300
301*     **** electron and geodesic data structures ****
302      call c_electron_init()
303      call c_geodesic_init()
304
305
306*     **** initialize QM/MM ****
307      call band_init_APC(rtdb)
308
309
310*     **** initialize HFX ****
311      call band_init_HFX(rtdb,nbrillioun,ispin,ne)
312
313
314*     **** initialize FixIon constraint ****
315      call ion_init_FixIon(rtdb)
316
317*     **** initialize linesearching ****
318      call linesearch_init()
319
320
321*                |**************************|
322******************   summary of input data  **********************
323*                |**************************|
324
325
326*     **** determine en ****
327      if (.not.control_spin_orbit()) then
328        icharge = 0.0d0
329        en(1)   = 0.0d0
330        en(2)   = 0.0d0
331        b = dble(3-cpsi_ispin())
332        do nb=1,brillioun_nbrillq()
333        weight = brillioun_weight(nb)
334        do ms=1,cpsi_ispin()
335          do i=1,ne(ms)
336            a = cpsi_occupation(nb,ms,i)
337            icharge = icharge - b*a*weight
338            en(ms)  = en(ms) + a*weight
339          end do
340        end do
341        end do
342        call K1dB_Vector_SumAll(2,en)
343        call K1dB_SumAll(icharge)
344      else
345        icharge          = -cpsi_ne(1)
346        en(1)            =  cpsi_ne(1)
347        en(cpsi_ispin()) =  cpsi_ne(cpsi_ispin())
348      end if
349
350
351      if (mprint) then
352         write(luout,1111) np
353         write(luout,1117) np_i,np_j,np_k
354         if (mapping.eq.1) write(luout,1112)
355         if (mapping.eq.2) write(luout,1113)
356         if (mapping.eq.3) write(luout,1118)
357
358         write(luout,1115)
359         write(luout,1121) control_boundry(),control_version()
360         if (.not.spin_orbit) then
361           if (cpsi_ispin().eq.1) write(luout,1130) "restricted"
362           if (cpsi_ispin().eq.2) write(luout,1130) "unrestricted"
363         else
364           write(luout,1130) "spin orbit"
365         end if
366
367         call v_bwexc_print(luout,control_gga())
368
369         call band_print_HFX(luout)
370         write(luout,1140)
371         do ia = 1,ion_nkatm()
372           write(luout,1150) ia,ion_atom(ia),
373     >                    cpsp_zv(ia),cpsp_lmax(ia)
374           comment = cpsp_comment(ia)
375           i = inp_strlen(comment)
376           write(6,1157) comment(1:i)
377           write(6,1158) cpsp_psp_type(ia)
378           write(luout,1152) cpsp_lmax(ia)
379           write(luout,1153) cpsp_locp(ia)
380c           write(luout,1154) cpsp_lmmax(ia)
381           write(luout,1154) cpsp_nprj(ia)
382           if (cpsp_semicore(ia))
383     >         write(luout,1155) cpsp_rcore(ia),cpsp_ncore(ia)
384           write(luout,1151) (cpsp_rc(i,ia),i=0,cpsp_lmax(ia))
385         end do
386
387         icharge = icharge + ion_TotalCharge()
388         write(luout,1159) icharge
389
390
391         write(luout,1160)
392         write(luout,1170) (ion_atom(K),ion_natm(K),K=1,ion_nkatm())
393         if (hprint) then
394         write(luout,1180)
395         do I=1,ion_nion()
396           if (ion_q_FixIon(I)) then
397           write(luout,1191) I,ion_aname(I),(ion_rion(K,I),K=1,3),
398     >                   ion_amass(I)/1822.89d0
399           else
400           write(luout,1190) I,ion_aname(I),(ion_rion(K,I),K=1,3),
401     >                   ion_amass(I)/1822.89d0
402           end if
403         end do
404         write(luout,1200) cx,cy,cz
405         write(luout,1210) gx,gy,gz
406         write(luout,1211) ion_nconstraints(),ion_ndof()
407         endif
408         write(luout,1220) en(1),en(cpsi_ispin()),
409     >                 ' (Fourier space)'
410         write(luout,1221) cpsi_ne(1),cpsi_ne(cpsi_ispin()),
411     >                 ' (Fourier space)'
412
413         write(luout,1230)
414         write(luout,1241) lattice_unita(1,1),
415     >                 lattice_unita(2,1),
416     >                 lattice_unita(3,1)
417         write(luout,1242) lattice_unita(1,2),
418     >                 lattice_unita(2,2),
419     >                 lattice_unita(3,2)
420         write(luout,1243) lattice_unita(1,3),
421     >                 lattice_unita(2,3),
422     >                 lattice_unita(3,3)
423         write(luout,1244) lattice_unitg(1,1),
424     >                 lattice_unitg(2,1),
425     >                 lattice_unitg(3,1)
426         write(luout,1245) lattice_unitg(1,2),
427     >                 lattice_unitg(2,2),
428     >                 lattice_unitg(3,2)
429         write(luout,1246) lattice_unitg(1,3),
430     >                 lattice_unitg(2,3),
431     >                 lattice_unitg(3,3)
432         call lattice_abc_abg(a,b,c,alpha,beta,gamma)
433         write(luout,1232) a,b,c,alpha,beta,gamma
434         write(luout,1231) lattice_omega()
435         write(luout,1260) ewald_rcut(),ewald_ncut()
436         write(luout,1261) ewald_mandelung()
437
438         write(luout,1255)
439         write(luout,1256) brillioun_nbrillioun(),
440     >                     brillioun_nbrillioun0()
441      end if
442
443c     **** print brillioun zone - extra logic for distributed kpoints ****
444      if (brillioun_print()) then
445         do i=1,brillioun_nbrillioun()
446            f0 = brillioun_weight_brdcst(i)
447            f1 = brillioun_ks_brdcst(1,i)
448            f2 = brillioun_ks_brdcst(2,i)
449            f3 = brillioun_ks_brdcst(3,i)
450            f4 = brillioun_k_brdcst(1,i)
451            f5 = brillioun_k_brdcst(2,i)
452            f6 = brillioun_k_brdcst(3,i)
453            if (mprint) write(luout,1257) f0,f1,f2,f3,f4,f5,f6
454         end do
455      else
456        if (mprint) write(luout,1258)
457      end if
458
459      if1 = Cram_nwave_all_brdcst(0)
460      if2 = Cram_nwave_brdcst(0)
461      if (mprint) then
462         write(luout,1249)
463         write(luout,1250) lattice_ecut(),ngrid(1),ngrid(2),ngrid(3),
464     >                     if1,if2
465      end if
466
467      if (brillioun_print()) then
468        do i=1,brillioun_nbrillioun()
469          if1 = Cram_nwave_all_brdcst(i)
470          if2 = Cram_nwave_brdcst(i)
471          if (mprint) then
472          write(luout,1251) i,lattice_wcut(),ngrid(1),ngrid(2),ngrid(3),
473     >                      if1,if2
474          end if
475        end do
476      else
477        if (mprint) write(luout,1252) lattice_wcut()
478      end if
479
480      if (mprint) then
481         if (control_smooth_cutoff())
482     >      write(luout,1253)
483     >      control_smooth_cutoff_values(1)*lattice_wcut(),
484     >      control_smooth_cutoff_values(2)
485
486         write(luout,1270)
487         write(luout,1280) control_time_step(),control_fake_mass()
488         write(luout,1290) control_tole(),control_tolc()
489         write(luout,1281) control_it_in()*control_it_out(),
490     >                 control_it_in(),control_it_out()
491
492         if ((minimizer.eq.5).or.(minimizer.eq.8)) then
493           write(6,1291)
494           write(6,1292)
495           write(luout,1295) control_ks_maxit_orb(),
496     >                       control_ks_maxit_orbs()
497           if (control_scf_algorithm().eq.0)
498     >       write(6,1293) "simple mixing"
499           if (control_scf_algorithm().eq.1)
500     >       write(6,1293) "Anderson potential mixing"
501           if (control_scf_algorithm().eq.2)
502     >       write(6,1293) "Johnson-Pulay mixing"
503           if (control_scf_algorithm().eq.3)
504     >       write(6,1293) "Anderson density mixing"
505           if (minimizer.eq.5) write(luout,1296) "potential"
506           if (minimizer.eq.8) write(luout,1296) "density"
507           write(6,1294) control_ks_alpha()
508           write(6,1301) control_kerker_g0()
509         end if
510         if (control_fractional()) then
511           write(6,1297)
512           if (control_fractional_smeartype().eq.0)
513     >       write(6,1298) "step function"
514           if (control_fractional_smeartype().eq.1)
515     >       write(6,1298) "Fermi-Dirac"
516           if (control_fractional_smeartype().eq.2)
517     >       write(6,1298) "Gaussian"
518           write(6,1299) control_fractional_kT(),
519     >                   control_fractional_temperature(),
520     >                   control_fractional_alpha()
521         end if
522         write(luout,1300)
523         call util_flush(luout)
524      end if
525
526
527
528*                |***************************|
529******************     call CG minimizer     **********************
530*                |***************************|
531      if (taskid.eq.MASTER) call current_second(cpu2)
532
533
534*     **** calculate energy ****
535      if (flag.eq.-1) then
536
537        EV= c_cgsd_noit_energy()
538      else
539
540        EV= c_cgsd_energy(newpsi)
541      end if
542
543*     **** calculate excited state orbitals ****
544      call ga_sync()
545      call c_cgsd_excited()
546
547*     **** extra energy output for QA test ****
548      if (lprint) write(luout,1600) EV
549
550*     **** calculate the spin contamination ****
551      if (flag.gt.-1) call cpsi_spin2(dipole(1))
552
553*     **** calculate the dipole ***
554      dipole(1) = 0.0d0
555      dipole(2) = 0.0d0
556      dipole(3) = 0.0d0
557
558*     ****  calculate gradient ***
559      if (flag.gt.0) then
560      nion = ion_nion()
561      value = BA_push_get(mt_dbl,(3*nion),
562     >                       'fion',fion(2),fion(1))
563      if (.not. value) call errquit('out of stack memory',0, MA_ERR)
564
565
566      call c_cgsd_energy_gradient(dbl_mb(fion(1)))
567      call dscal(3*nion,(-1.0d0),dbl_mb(fion(1)),1)
568      end if
569
570*     **** calculate the stress tensor ****
571      call dcopy(9,0.0d0,0,stress,1)
572      call dcopy(6,0.0d0,0,lstress,1)
573      if (flag.eq.3) then
574
575         call cpsp_stress_init()
576         call cpsp_stress_readall()
577         call c_cgsd_energy_stress(stress,lstress)
578         call cpsp_stress_end()
579      end if
580
581
582*     *************************************************************
583*     **** output energy, dipole, and gradient to rtdb for use ****
584*     **** by task_energy and task_gradient                    ****
585*     *************************************************************
586      value = .true.
587      if (flag.gt.-1) then
588      value = btdb_put(rtdb,'band:energy',mt_dbl,1,EV)
589      value = value.and.
590     >        btdb_put(rtdb,'band:dipole',mt_dbl,
591     >                 3,dipole)
592      end if
593      if (flag.gt.0) then
594        value = value.and.btdb_put(rtdb,'band:gradient',mt_dbl,
595     >                             3*nion,dbl_mb(fion(1)))
596        value = value.and.BA_pop_stack(fion(2))
597      end if
598      if (flag.eq.3) then
599        value = value.and.
600     >        btdb_put(rtdb,'band:stress',mt_dbl,
601     >                 9,stress)
602        value = value.and.
603     >        btdb_put(rtdb,'band:lstress',mt_dbl,
604     >                 6,lstress)
605      end if
606      if (.not. value)
607     > call errquit('band_minimizer: error writing rtdb',0, RTDB_ERR)
608
609
610
611      if (taskid.eq.MASTER) call current_second(cpu3)
612
613*                |***************************|
614******************         Epilogue          **********************
615*                |***************************|
616
617*     **** calculate interpolated band structure plot ***
618      call band_interpolate_structure(rtdb)
619cc      call band_kp_structure(rtdb)
620
621*     **** calculate Mulliken Populations ***
622      if (control_Mulliken()) then
623         call cpsi_Mulliken(rtdb)
624      end if
625
626
627*     **** write geometry to rtdb ****
628      call ion_write(rtdb)
629
630*     **** write wavefunctions to file and finalize psi ****
631      if (flag.eq.-1) then
632        value = cpsi_finalize(.false.)
633      else
634        value = cpsi_finalize(.true.)
635      end if
636
637*     **** deallocate heap memory ****
638      call c_electron_finalize()
639      call c_geodesic_finalize()
640      call ewald_end()
641      call cstrfac_end()
642      call c_coulomb_end()
643      call cke_end()
644      call cpsp_end()
645      call Cram_end()
646      call c_G_end()
647      call band_end_HFX()
648      call brillioun_end()
649      call band_end_APC()
650      call ion_end()
651      call ion_end_FixIon()
652      call cpsi_data_end()
653      call C3dB_pfft_end()
654      call C3dB_end(1)
655      IF ((control_gga().ge.10).and.(control_gga().le.200)) THEN
656      call mask_end()
657      call G_end()
658      call D3dB_end(1)
659      end if
660
661*                |***************************|
662****************** report consumed cputime   **********************
663*                |***************************|
664      if (lprint) then
665         CALL current_second(cpu4)
666
667         T1=CPU2-CPU1
668         T2=CPU3-CPU2
669         T3=CPU4-CPU3
670         T4=CPU4-CPU1
671         AV=T2/dble(c_electron_count())
672         write(luout,1801)
673         write(luout,1802)
674         write(luout,1803) T1
675         write(luout,1804) T2
676         write(luout,1805) T3
677         write(luout,1806) T4
678         write(luout,1807) AV,c_electron_count(),linesearch_count()
679
680         call nwpw_timing_print_final(mprint,c_electron_count())
681         write(luout,*)
682         CALL nwpw_MESSAGE(4)
683      end if
684
685
686      call Parallel3d_Finalize()
687      call Parallel_Finalize()
688      band_minimizer = value
689      return
690
691
692*:::::::::::::::::::::::::::  format  :::::::::::::::::::::::::::::::::
693 1000 FORMAT(10X,
694     > '**********************************************************')
695 1010 FORMAT(10X,
696     > '*                                                        *')
697 1020 FORMAT(10X,
698     > '*                   NWPW BAND Calculation                *')
699 1030 FORMAT(10X,
700     > '*  [(bundled Grassmann/Stiefel manifold implementation)] *')
701 1035 FORMAT(10x,
702     > '*         [ NorthWest Chemistry implementation ]         *')
703 1040 FORMAT(10X,
704     > '*                version #1.10   01/31/03                *')
705 1041 FORMAT(10X,
706     > '*  A pseudopotential plane-wave band structure program   *')
707 1042 FORMAT(10X,
708     > '*  with Brillouin zone sampling for optimizing crystals, *')
709 1043 FORMAT(10X,
710     > '*  slabs, and polymers.  Developed by Eric J. Bylaska,   *')
711 1044 FORMAT(10X,
712     > '*  Edoardo Apra, and Patrick Nichols.                    *')
713 1100 FORMAT(//)
714 1110 FORMAT(10X,'================ input data ========================')
715 1111 FORMAT(/' number of processors used:',I16)
716 1112 FORMAT( ' parallel mapping         :         1d slab')
717 1113 FORMAT( ' parallel mapping         :     2d  hilbert')
718 1115 FORMAT(/' options:')
719 1117 FORMAT( ' processor grid           :',I4,' x',I4,' x',I4)
720 1118 FORMAT( ' parallel mapping         :       2d hcurve')
721 1120 FORMAT(5X,' ionic motion         = ',A)
722 1121 FORMAT(5X,' boundary conditions  = ',A,'(version', I1,')')
723 1130 FORMAT(5X,' electron spin        = ',A)
724 1131 FORMAT(5X,' exchange-correlation = ',A)
725 1140 FORMAT(/' elements involved in the cluster:')
726 1150 FORMAT(5X,I2,': ',A4,'  core charge:',F4.1,'  lmax=',I1)
727 1151 FORMAT(5X,'        cutoff =',4F8.3)
728 1152 FORMAT(12X,' highest angular component      : ',i3)
729 1153 FORMAT(12X,' local potential used           : ',i3)
730 1154 FORMAT(12X,' number of non-local projections: ',i3)
731 1155 FORMAT(12X,' semicore corrections included  : ',
732     >       F6.3,' (radius) ',F6.3,' (charge)')
733 1156 FORMAT(12X,' aperiodic cutoff radius        : ',F6.3)
734 1157 FORMAT(12X,' comment    : ',A)
735 1158 FORMAT(12X,' pseudpotential type            : ',i3)
736
737 1159 FORMAT(/' total charge:',F8.3)
738 1160 FORMAT(/' atomic composition:')
739 1170 FORMAT(7(5X,A4,':',I3))
740 1180 FORMAT(/' initial position of ions:')
741 1190 FORMAT(5X, I4, A5  ,' (',3F11.5,' ) - atomic mass= ',F7.3,' ')
742 1191 FORMAT(5X, I4, A5, ' (',3F11.5,
743     >       ' ) - atomic mass= ',F7.3,' - fixed')
744 1200 FORMAT(5X,'   G.C.  ',' (',3F11.5,' )')
745 1210 FORMAT(5X,'   C.O.M.',' (',3F11.5,' )')
746 1211 FORMAT(5X,'   number of constraints = ', I6,' ( DOF = ',I6,' )' )
747 1220 FORMAT(/' number of electrons: spin up=',F8.2,
748     >                          '  spin down=',F8.2,A)
749 1221 FORMAT( ' number of orbitals:  spin up=',I8,
750     >                          '  spin down=',I8,A)
751 1230 FORMAT(/' supercell:')
752 1231 FORMAT(5x,' volume : ',F10.1)
753 1232 FORMAT(/5x,' lattice:    a=',f8.3,'    b=',f8.3,'     c=',f8.3,
754     >       /5x,'         alpha=',f8.3,' beta=',f8.3,' gamma=',f8.3)
755 1241 FORMAT(5x,' lattice:    a1=<',3f8.3,' >')
756 1242 FORMAT(5x,'             a2=<',3f8.3,' >')
757 1243 FORMAT(5x,'             a3=<',3f8.3,' >')
758 1244 FORMAT(5x,' reciprocal: b1=<',3f8.3,' >')
759 1245 FORMAT(5x,'             b2=<',3f8.3,' >')
760 1246 FORMAT(5x,'             b3=<',3f8.3,' >')
761
762 1249 FORMAT(/' computational grids:')
763 1250 FORMAT(5X,' density     cutoff=',F7.3,'  fft=',I4,'x',I4,'x',I4,
764     &       '( ',I8,' waves ',I8,' per task)')
765 1251 FORMAT(5X,' wavefnc ',I3,' cutoff=',F7.3,
766     &        '  fft=',I4,'x',I4,'x',I4,
767     &       '( ',I8,' waves ',I8,' per task)')
768 1252 FORMAT(5x,' wavefnc     cutoff=',F7.3,
769     >       ' wavefunction grids not printed - ',
770     >       'number of k-points is very large')
771 1253 FORMAT(5X,' smooth wavefnc cutoff:  A=',F7.3,' sigma = ',F7.3)
772
773
774 1255 FORMAT(/' brillouin zone:')
775 1256 FORMAT(5x,' number of zone points:',I6,
776     >          ' (reduced from ',I0,' zone points without symmetry)')
777 1257 FORMAT(5x,' weight=',f8.3,'  ks=<',3f8.3,' >, k=<',3f8.3,'>')
778 1258 FORMAT(5x,' number of k-points is very large')
779
780 1260 FORMAT(5X,' Ewald summation: cut radius=',F8.2,'  and',I3)
781 1261 FORMAT(5X,'                   madelung=',f14.8)
782
783 1270 FORMAT(/' technical parameters:')
784 1280 FORMAT(5X, ' time step=',F10.2,5X,'fictitious mass=',F10.1)
785 1281 FORMAT(5X, ' maximum iterations =',I10,
786     >           ' ( ',I4,' inner ',I6,' outer )')
787 1290 FORMAT(5X, ' tolerance=',E8.3,' (energy)',E12.3,
788     &        ' (density)')
789 1291 FORMAT(/' Kohn-Sham scf parameters:')
790 1292 FORMAT(5X, ' Kohn-Sham algorithm  = conjugate gradient')
791 1293 FORMAT(5X, ' scf algorithm        = ',A)
792 1294 FORMAT(5X, ' scf mixing parameter =',F7.4)
793 1295 FORMAT(5X, ' Kohn-Sham iterations = ',I3,
794     >           ' (',I3,' outer)')
795 1296 FORMAT(5X, ' scf mixing type      = ',A)
796 1297 FORMAT(/' fractional smearing parameters:')
797 1298 FORMAT(5X, ' smearing algorithm   = ',A)
798 1299 FORMAT(5X, ' smearing parameter   = ',E9.3,' (',F7.1,' K)'/,
799     >       5X, ' mixing parmameter    = ',F7.4)
800 1300 FORMAT(//)
801 1301 FORMAT(5X, ' Kerker damping       =',F7.4)
802 1305 FORMAT(10X,'================ iteration =========================')
803 1310 FORMAT(I8,E20.10,3E15.5)
804 1320 FORMAT(' number of electrons: spin up=',F11.5,'  down=',F11.5,A)
805 1330 FORMAT(/' comparison between hamiltonian and lambda matrix')
806 1340 FORMAT(I3,2I3,' H=',E16.7,', L=',E16.7,', H-L=',E16.7)
807 1350 FORMAT(/' orthonormality')
808 1360 FORMAT(I3,2I3,E18.7)
809 1370 FORMAT(I3)
810 1380 FORMAT(' ''',a,'''',I4)
811 1390 FORMAT(I3)
812 1400 FORMAT(I3,3E18.8/3X,3E18.8)
813 1410 FORMAT(10X,'=============  summary of results  =================')
814 1420 FORMAT( ' final position of ions:')
815 1430 FORMAT(/' total     energy    :',E19.10,' (',E15.5,'/ion)')
816 1440 FORMAT( ' total orbital energy:',E19.10,' (',E15.5,'/electron)')
817 1450 FORMAT( ' hartree   energy    :',E19.10,' (',E15.5,'/electron)')
818 1460 FORMAT( ' exc-corr  energy    :',E19.10,' (',E15.5,'/electron)')
819 1470 FORMAT( ' ion-ion   energy    :',E19.10,' (',E15.5,'/ion)')
820 1480 FORMAT(/' K.S. kinetic energy :',E19.10,' (',E15.5,'/electron)')
821 1490 FORMAT( ' K.S. V_l  energy    :',E19.10,' (',E15.5,'/electron)')
822 1495 FORMAT( ' K.S. V_nl energy    :',E19.10,' (',E15.5,'/electron)')
823 1496 FORMAT( ' K.S. V_Hart energy  :',E19.10,' (',E15.5,'/electron)')
824 1497 FORMAT( ' K.S. V_xc energy    :',E19.10,' (',E15.5,'/electron)')
825 1498 FORMAT( ' Virial Coefficient  :',E19.10)
826 1500 FORMAT(/' orbital energies:')
827 1510 FORMAT(2(E18.7,' (',F8.3,'eV)'))
828 1600 FORMAT(/' Total BAND energy   :',E19.10)
829
830 1801 FORMAT(//'== Timing ==')
831 1802 FORMAT(/'cputime in seconds')
832 1803 FORMAT( '  prologue    : ',E14.6)
833 1804 FORMAT( '  main loop   : ',E14.6)
834 1805 FORMAT( '  epilogue    : ',E14.6)
835 1806 FORMAT( '  total       : ',E14.6)
836 1807 FORMAT( '  cputime/step: ',E14.6,
837     >        '       (',I8,' evalulations,', I8,' linesearches)')
838 1808 FORMAT(A,E14.6,E14.6)
839 1809 FORMAT(//A,2A14)
840
841 9010 FORMAT(//' >> job terminated due to code =',I3,' <<')
842
843 9000 if (taskid.eq.MASTER) write(6,9010) ierr
844      call Parallel_Finalize()
845
846      band_minimizer = value
847      return
848      END
849
850