******************************************************************************** * * * 2DHF version 1-2003 * * Copyright (C) 1996 Jacek Kobus, Leif Laaksonen, Dage Sundholm * * * * This software may be used and distributed according to the terms * * of the GNU General Public License, see README and COPYING. * * * ******************************************************************************** c ### scf #### c c Controls SCF process. c subroutine scf (cw_sor,cw_orb,cw_coul,cw_exch,cw_suppl,cw_sctch) implicit integer*4 (i-n) implicit real*8 (a-h,o-z) logical stop_x2dhf integer*4 cw_sor character*1 str(80) c_mpi include 'mpif.h' c_mpi include 'commons8.inc' include 'commons_mpi.inc' dimension cw_orb(*),cw_coul(*),cw_exch(*), & cw_suppl(*),cw_sctch(*) dimension inde(60),demaxt(60) equivalence(datetime,str(1)) parameter (msglen=1000000000) call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr ) c prepare scf process c normalize/orthogonalize orbitals, calculate Lagrange multipliers c prepare data for boundary value evaluation c if (myid.eq.0) then call prepscf (cw_sor,cw_orb,cw_coul,cw_exch,cw_suppl,cw_sctch) endif c broadcast orbitals and potentials call MPI_BCAST(cw_orb,length1,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) c print *,'about to broadcast cw_coul' call MPI_BCAST(cw_coul,length2,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) c print *,'about to broadcast cw_exch' if (length3.ge.msglen) then length3a=length3/2 length3b=length3-length3a call MPI_BCAST(cw_exch,length3a,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) call MPI_BCAST(cw_exch(length3a+1),length3b, & MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) else call MPI_BCAST(cw_exch,length3,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) endif stop_x2dhf=.false. ic=0 thren=10.d0**(-ienterm) thrno=10.d0**(-inoterm) c the scf procedure starts here (set some parameters and go ...) ddmax=0.d0 do i=1,norb inde(i)=i itouch(i)=0 demax(i)=100.d0 engi(i)=engo(i) enddo dnomax=100.d0 tscf=0.d0 trelax=0.d0 tortho=0.d0 trayl =0.d0 tlagra=0.d0 tmomen=0.d0 c if extension has not beeen read in recalculated asymptotic values c for potentials if (idump.ne.0.and.(iform.eq.0.or.iform.eq.3)) iasympt=norb call ttime (time1) tscf=time1 iscf=0 time1=0.d0 time2=0.d0 if (ihomon.eq.2) then call homo(0,cw_orb) endif c initialize totenp if (myid.eq.0) then c total energy is only calculated by rank 0 process if (iform.eq.0.or.iform.eq.2) then iorb=0 call totenp (iorb,cw_orb,cw_coul,cw_exch, & cw_suppl(i4b( 4)),cw_suppl(i4b( 5)),cw_suppl(i4b(13)), & cw_suppl(i4b(14)), & cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),cw_sctch(i5b( 3)), & cw_sctch(i5b( 4))) endif endif c ----- begin of scf loop ----- if (myid.eq.0) then if (maxscf.le.0) write(iout6,*) & '... skipping scf iterations ...' endif do iscf=1,maxscf if (myid.eq.0) then if (iscf.eq.1) write(iout6,15000) endif c loop over orbitals (in reverse order) c relaxation begins from the innermost orbital, i.e. the last c one on the input data list ibeg=1 iend=norb c prepare retrieval of exchange potentials c this could be done once in front of the loop call prepwexch c ----- begin of orbital loop ----- do i=ibeg,iend iorb=norb-i+1 itouch(iorb)=1 if (ifix(iorb).ne.0) goto 100 call ttime (time1) if (myid.eq.0) then if (iform.eq.0.or.iform.eq.2) then c retrieve from disk exchange potentials involving iorb call rfdexch(iorb,cw_exch) if (iasympt.ne.0) then call asymptp(iorb,cw_coul,cw_exch) iasympt=iasympt-1 endif endif endif call MPI_BCAST(i3b,1830,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) call MPI_BARRIER(MPI_COMM_WORLD,ierr ) if (length3.ge.msglen) then length3a=length3/2 length3b=length3-length3a call MPI_BCAST(cw_exch,length3a,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) call MPI_BCAST(cw_exch(length3a+1),length3b, & MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) else c call dcopy(length3,cw_exch,1,cw_exch_mpi,1) c call MPI_BCAST(cw_exch,length3,MPI_DOUBLE_PRECISION,0, c & MPI_COMM_WORLD,ierr) c call MPI_BCAST(cw_exch_mpi,length3,MPI_DOUBLE_PRECISION, call MPI_BCAST(cw_exch,length3,MPI_DOUBLE_PRECISION, & 0,MPI_COMM_WORLD,ierr) endif call MPI_BARRIER( MPI_COMM_WORLD,ierr ) c write(*,'("before: ",2i4,(5d15.6))') myid,iorb, c & (cw_exch(i1),i1=100,length3,nni*mxnmu) call relax (iorb,cw_sor,cw_orb,cw_coul,cw_exch, & cw_suppl(i4b( 1)),cw_suppl(i4b( 2)),cw_suppl(i4b( 3)), & cw_suppl(i4b( 4)),cw_suppl(i4b( 5)),cw_suppl(i4b( 6)), & cw_suppl(i4b( 7)),cw_suppl(i4b( 8)),cw_suppl(i4b( 9)), & cw_suppl(i4b(10)),cw_suppl(i4b(11)),cw_suppl(i4b(12)), & cw_suppl(i4b(14)), & cw_sctch(i5b( 1)),cw_sctch(i5b( 2)),cw_sctch(i5b( 3)), & cw_sctch(i5b( 4)),cw_sctch(i5b( 5)),cw_sctch(i5b( 6))) c write(*,'("after: ",2i4,(5d15.6))') myid,iorb, c & (cw_exch(i1),i1=100,length3,nni*mxnmu) call ttime (time2) trelax=trelax+(time2-time1) if (myid.eq.0) then time1=time2 call norm (iorb,cw_orb, & cw_suppl(i4b(9)),cw_suppl(i4b(14)),cw_sctch(i5b(1))) c enforce C_i symmetry of the orbital if (ihomon.eq.2) call homo(iorb,cw_orb) iprt=0 if (norb.gt.1) then call ortho (iorb,cw_orb, & cw_suppl(i4b(9)),cw_suppl(i4b(14)), & cw_sctch(i5b(1)),iprt) endif call ttime(time2) tortho=tortho+(time2-time1) c calculate diagonal energy value call ttime(time1) call rayl (iorb,cw_orb,cw_coul,cw_exch, & cw_suppl(i4b( 4)),cw_suppl(i4b( 5)), & cw_suppl(i4b(13)),cw_suppl(i4b(14)), & cw_sctch(i5b( 1)),cw_sctch(i5b( 2)), & cw_sctch(i5b( 3)),cw_sctch(i5b( 4))) demaxt(iorb)=eng(iorb)-engi(iorb) engi(iorb)=eng(iorb) c calculate off-diagonal Lagrange multipliers call lagra (iorb,cw_orb,cw_coul,cw_exch, & cw_suppl(i4b(14)),cw_sctch(i5b(1)),cw_sctch(i5b(2))) call ttime(time2) trayl =trayl + (time2-time1) endif c demaxt call MPI_BCAST(demaxt,60,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) c eng,engi call MPI_BCAST(eng,120,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) c common/noeng/engo(3600),h(1,1),inon,inpv(10),icou(60,60) call MPI_BCAST(engo,3600,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) call MPI_BCAST(area,60,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) c broadcast updated orbitals call MPI_BCAST(cw_orb,length1,MPI_DOUBLE_PRECISION,0, & MPI_COMM_WORLD,ierr) if (myid.eq.0) then c save updated exchange potentials if (iform.eq.0.or.iform.eq.1) then call wtdexch(iorb,cw_exch) endif if (iprtlev.eq.1) then write(iout6,15100) iscf,iorn(iorb),bond(iorb), & gut(iorb),demaxt(iorb),eng(iorb),area(iorb) endif endif if (myid.eq.0) then if (nobckup.gt.0) then modv=mod(iscf,iabs(nobckup)) if (modv.eq.0) then if (iform.eq.0.or.iform.eq.2) then c calculate 'intermidiate' total energy call totenp (iorb,cw_orb,cw_coul,cw_exch, & cw_suppl(i4b( 4)),cw_suppl(i4b( 5)), & cw_suppl(i4b(13)),cw_suppl(i4b(14)), & cw_sctch(i5b( 1)),cw_sctch(i5b( 2)), & cw_sctch(i5b( 3)),cw_sctch(i5b( 4))) endif endif endif endif c Sometimes a lengthy calculation needs to be stoped earlier. c kill will do the job but you risk that the data being written c will be lost. To avoid this and force the program to end gracefully c make c touch stop_x2dhf c If this file is detected the program ends. inquire (file ='stop_x2dhf',exist=stop_x2dhf) if (stop_x2dhf) goto 110 100 continue enddo c ----- end of orbital loop ----- 110 continue c if (myid.eq.0) then demax(1)=0.d0 inde(1)=1 dnomax=0.d0 if (norb.gt.0) then do nei=1,norb if(abs(demaxt(nei)).ge.abs(demax(1))) then demax(1)=demaxt(nei) inde(1)=nei endif enddo do nei=1,norb if (abs(1.d0-area(nei)).ge.dnomax) then dnomax=abs(1.d0-area(nei)) endif enddo endif c inde containes indexes of orbitals in the order of their state c of convergence; inde(1) points to the worst converged orbital if (ddmax.eq.0.d0) ddmax=demax(1) if (myid.eq.0) then if (iprtlev.eq.2) then write(iout6,15110) iscf,iorn(inde(1)),bond(inde(1)), & gut(inde(1)),demaxt(iorb),eng(iorb),area(iorb) endif endif c check if the scf process is to be continued if (abs(demax(1)).gt.diver) then write(*,*) ' ' write(*,*) '... maximal orbital energy error is too ', & 'large ...' goto 9999 elseif (iscf+1.gt.maxscf) then write(*,*) ' ' write(*,*) '... scf iteration limit reached ... ' goto 9999 elseif (abs(demax(1)).lt.thren) then write(*,*) ' ' write(*,*) '... orbital energy threshold reached ...' goto 9999 elseif (exlorb.eq.0.and.dnomax.lt.thrno) then write(*,*) ' ' write(*,*) '... orbital normalization threshold', & ' reached ...' goto 9999 elseif (stop_x2dhf) then write(*,*) ' ' write(*,*) '... stop_x2dhf detected ... ' write(*,*) '... aborting gracefully ...' write(*,*) goto 9999 endif c write functions to disk file if nobckup > 0 if (myid.eq.0) then if (nobckup.gt.0) then modv=mod(iscf,nobckup) if (modv.eq.0) then call wtdisk(cw_orb,cw_coul,cw_exch) if (iprtlev.eq.3) write(iout6,15110) & iscf,iorn(inde(1)),bond(inde(1)),gut(inde(1)), & demaxt(iorb),eng(iorb),area(iorb) if (iform.eq.1.or.iform.eq.3) then c calculate 'intermidiate' total energy call toten (cw_orb,cw_coul,cw_exch, & cw_suppl(i4b( 4)),cw_suppl(i4b( 5)), & cw_suppl(i4b(13)),cw_suppl(i4b(14)), & cw_sctch(i5b( 1)),cw_sctch(i5b( 2)), & cw_sctch(i5b( 3)),cw_sctch(i5b( 4))) endif write(iout6,6100) etot endif endif endif c let all processes calculate multipole moments c recalculate multipole moments expansion coefficients if (abs(ddmax/demax(1)).gt.facmul.or. & abs(demax(1)/ddmax).gt.facmul) then call momen(cw_orb,cw_suppl,cw_sctch) ddmax=demax(1) c recalculate asymptotic values of potentials if (iform.eq.1.or.iform.eq.3) then call asympot (cw_coul,cw_exch) else iasympt=norb endif do iorb=1,norb itouch(iorb)=0 enddo endif c to run a job in a bach mode (or in the background) with c the output redirected into a file and to monitor every single c iteration if (myid.eq.0) then call flush(iout6) endif enddo c ----- end of scf loop ----- 09999 iscf=iscf-1 if (myid.eq.0) then call summary (iscf,cw_orb,cw_coul,cw_exch,cw_suppl,cw_sctch) call ttime (time1) tscf=time1-tscf write(*,5010) tscf call getdatetime(datetime) call clrstring(imax,datetime) write(*,'(//,1x,80a1)') (str(i),i=1,imax) endif write(*,1305) myid,trelax 1300 format(/,80a1) 1305 format(/,'pe =',i2,' cpu in relax ',f12.2) 5010 format(4x,'total time in scf routine ', f12.2) 6100 format(1x,'total energy: ',d20.13) 15000 format(/' iter. orbital',7x,'energy diff.',6x,'orbital energy', & 12x,'orbital norm') 15100 format(i5,i4,1x,a8,a1,3x,d13.5,2d25.14) 15110 format(i5,i4,1x,a8,a1,3x,d13.5,2d25.14) return end