Opium v4.0

rappegroup · Mar 26, 2021 · c81be50 · c81be50
1 parent 9393e72
commit c81be50
Show file tree

Hide file tree

Showing 375 changed files with 23,788 additions and 514,170 deletions.
diff --git a/install.command b/install.command
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd "$(dirname "$0")"
+./configure
+make
+#sudo cp opium /usr/local/bin/opium
+echo "*****************************************************************************"
+echo "Installation finished! Run opium using command line in a new terminal window!"
+echo "The usual location for the executable is /usr/local/bin"
+echo "*****************************************************************************"
+#make distcheck
+# make install
diff --git a/src/atom/average.F b/src/atom/average.F
@@ -48,6 +48,8 @@ subroutine average(ixc,exccut,iexit)
       common /iavgtype/ iavgtype
 
       dimension rvxct(npdm),rexct(npdm),v(npdm)
+      dimension rvxt(npdm),rext(npdm)
+      dimension rvct(npdm),rect(npdm)
       dimension ibp(n0),indall(n0),rsatom(npdm),rv(npdm)
       dimension fe_old(n0),fe(n0), rnlt(npdm,n0)
       dimension rvscr(npdm),rvh(npdm),phom(npdm)
@@ -322,7 +324,8 @@ subroutine average(ixc,exccut,iexit)
          enddo
 
 c     3a: Find new xc and hartree potentials since rho has changed
-         call excorr (np,ixc,exccut,rsatom,rvxct,rexct,irelxc)
+         call excorr (np,ixc,exccut,rsatom,rvxct,rexct,
+     $        rvxt,rext,rvct,rect,irelxc)
          call hrtree (np,h,r,rsatom,rvh)
 
 c     rvscr is the screening potential
@@ -556,7 +559,8 @@ subroutine average(ixc,exccut,iexit)
       call radin(r,rsatom,0,maxim,h,pow)
       pow = 2.0
       call radin(r,rscore,0,maxim,h,pow)
-      call excorr (maxim,ixc,exccut,rsatom,rvxct,rexct,irelxc)
+      call excorr (maxim,ixc,exccut,rsatom,rvxct,rexct,
+     $     rvxt,rext,rvct,rect,irelxc)
       call hrtree (maxim,h,r,rsatom,rvh)
       do  i = maxim+1,np
         rvh(i)  = z + z - xion - xion

diff --git a/src/atom/dftsolve.F b/src/atom/dftsolve.F
@@ -16,6 +16,7 @@
 c Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 c
 c
+c change some of the 911 case to be less than 50 instead of less than 10
       subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
      $     irelxc,iprint)
 
@@ -58,26 +59,31 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
       common /aval/ rcall(n0),rvap(n0),rnorm(n0),ibd(n0),etot
 
       common /niter/ niter
+      common /ahf/ s(n0),az(n0),acc(n0),dpm(n0)
 c     -------------------------------------------------------------------------
 c     Internal (Fortran only) common blocks                
 c     -------------------------------------------------------------------------
       common /nlpot1/ Flstar(npdm),phipsref(npdm,n0),phiguess(npdm,n0)
-      common /partpot/ rvh(npdm),rvxc(npdm),rexc(npdm)
+      common /partpot/ rvh(npdm),rvxc(npdm),rexc(npdm),
+     $     rvx(npdm),rex(npdm),rvc(npdm),rec(npdm),rvdft(npdm)
       common /nlcore/ rvloc(npdm)
       common /ipos/ ipos(n0),itermcount
       common /vdir/ vdir(npdm),pold(npdm)
       common /rold/ rnlold(npdm,n0),rpold(npdm,n0)
+      common /ameth/ meth(n0),imeth,nfails
 c     -------------------------------------------------------------------------
 
 c     *************************************************************************
 c     local variables
 c     *************************************************************************
       dimension rvn(npdm),rvf(npdm),eold(n0),p(npdm)
-      dimension rsatom(npdm),f(npdm),f2(npdm),f3(npdm)
+      dimension rsatom(npdm),f(npdm),f2(npdm),f3(npdm),f4(npdm)
       dimension pder(npdm),ader(npdm)
       dimension dl(n0), rsatom2(npdm)
       dimension rsold(npdm),wnlfract(npdm)
       dimension wavea(npdm),icarray(n0)
+      dimension emattest(n0,n0)
+      dimension lcount(10)
 c     rvpostloc(npdm),rvpreloc(npdm)
 c     *************************************************************************
 c     Section 1:  set up tolerance params and zero out arrays.
@@ -176,6 +182,7 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
         rvcoul(i) = rvps(i,1) - rvcore(i,1)
       enddo
 
+
       if (inl.ne.0) then
         write(7,*)"Performing non-local pseudopotential calculation"
       else
@@ -248,6 +255,11 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
          npratt = npratt + 1
          if (npratt.gt.ipratt) npratt = 0
 
+c         do i = 1,np
+c            write(80,*) rsatom(i)
+c         enddo
+
+
          edmax = 0
          difmax = 0
          maxim = 0
@@ -262,6 +274,7 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
 c     ***********************************************************************
 
          do m = 1,norb
+
             if (m.le.ncore) then
                ibound=1
             else
@@ -276,7 +289,7 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
 
                   if (inl.ne.0) then
 c     if non-local
-                     call applyaug(rvloc,rvps(1,iloc))                        
+                     call applyaug(rvloc,rvps(1,iloc))
                      do i=1,np
                         Flstar(i) = (rvps(i,m) - rvloc(i))
      $                       *phipsref(i,m)/r(i)
@@ -285,7 +298,9 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
                      call schsl (m,no(m),lo(m),en(m),nmax(m),
      $                    rvloc,p,igh,ibound)
 
-                     if (itermcount.gt.10) goto 911
+                     if (itermcount.gt.50) then
+                         goto 911
+                     endif 
 c     911 = terminal error - exit scpot
 
                      call applyaug(rvloc,rvcore(1,iloc))
@@ -296,11 +311,13 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
 
                   else
 c     if AE or local 
-
                      call schsl(m,no(m),lo(m),en(m),nmax(m),
      $                    rvps(1,m),p,igh,ibound)
 
-                     if (itermcount.gt.10) goto 911
+                     if (itermcount.gt.50) then
+                         write(7,*) "line 445"
+                         goto 911
+                     endif
 c     911 = terminal error - exit scpot
 
                   endif
@@ -327,7 +344,9 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
                   call schslrel(m,rvps(1,m),rnla(1,m),
      $                 rnlb(1,m),z,rvcoul,en,nmax(m),ibound)
 
-                  if (itermcount.gt.10) goto 911
+                  if (itermcount.gt.50) then
+                      goto 911
+                  endif
 c     911 = terminal error - exit scpot
 
                   do i=1,np
@@ -383,21 +402,24 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
             do i=1,maxim
                rsatom2(i)=rsatom(i)+rscore(i)
             enddo
-            call excorr(maxim,ixc,exccut,rsatom2,rvxc,rexc,irelxc)
+            call excorr(maxim,ixc,exccut,rsatom2,
+     $           rvxc,rexc,rvx,rex,rvc,rec,irelxc)
 
             do i = 1,maxim
                f(i) = 0.5d0 * rvh(i)*rsatom(i)/r(i)
      $              + rexc(i)*rsatom2(i)/r(i)
      $              - (rvxc(i)*rsatom(i)+rvh(i)*rsatom(i))/r(i)
             enddo
          else
-            call excorr(maxim,ixc,exccut,rsatom,rvxc,rexc,irelxc)
+            call excorr(maxim,ixc,exccut,rsatom,
+     $           rvxc,rexc,rvx,rex,rvc,rec,irelxc)
 
             do i = 1,maxim
                f(i) = (0.5d0 * rvh(i) + rexc(i) - (rvxc(i)+rvh(i)) ) 
      $              * rsatom(i)/r(i)
                f2(i) = 0.5d0 * rvh(i)* rsatom(i)/r(i)
                f3(i) = rexc(i)* rsatom(i)/r(i)
+               f4(i) = rvxc(i)* rsatom(i)/r(i)
             enddo
          endif
          ehxc = 2
@@ -406,11 +428,14 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
          call radin(r,f2,0,maxim,h,ehar)
          exc = 2
          call radin(r,f3,0,maxim,h,exc)
+         evxc = 2
+         call radin(r,f4,0,maxim,h,evxc)
          ebs = 0
          do m = 1,norb
             ebs = ebs + wnl(m) * en(m)
          enddo
          etot = ebs + ehxc
+         etot = ebs - ehar + exc - evxc
 
 c     ***********************************************************************
 c     Section 6:  Mix potentials
@@ -419,7 +444,7 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
          dvmax = 0
          xn2 = (zeff - xion) * 2
          do i = 1,maxim
-            coul = rvh(i) + rvxc(i)
+            coul = rvh(i) + rvxc(i)  
             dv = abs(coul - rvcoul(i))
             if (dv.gt.dvmax) igdd=i
             dvmax = max(dv,dvmax)
@@ -443,7 +468,9 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
      $              + (1.0 - xmix) * rvps(i,iorb)
             enddo
          enddo
-         if (niter.eq.maxit) goto 911
+         if (niter.eq.maxit) then
+             goto 911
+          endif
 
 c     911 = terminal error - exit scpot
 
@@ -512,7 +539,7 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
                   write (7,505) nlm(m),en(m)               
 
                   if (en(m).ne.en(m)) then
-                     goto 911
+                     goto 911 
 c     911 = terminal error - exit scpot
                   endif
                endif
@@ -521,6 +548,11 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
 
       enddo
 
+c      close(77)
+c      close(78)
+c      close(79)
+c      close(80)
+
  505  format(1x,'!WARNING! positive eigenvalue predicted',
      $     ' in scpot  nlm=',i4,5x,'new e =',e13.7)
 
@@ -814,7 +846,7 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
       endif
 
       write(7,701) etot,ebs,ehxc
-      write(7,702) ehar,exc
+      write(7,702) ehar,exc,evxc
 
       write (7,*)
       write (7,9222)
@@ -876,7 +908,7 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
  500  format(1x,i4,1x,'iterations: (units are Ry and bohr)')
  700  format(1x,'Converged in 1 iteration (probably reference state)')
  701  format(1x,'Energy: ',f16.8,2x,'Ebs: ',f16.8,2x,'Ehxc: ',f16.8)
- 702  format(1x,'Eh    : ',f16.8,2x,'Exc: ',f16.8)
+ 702  format(1x,'Eh    : ',f16.8,2x,'Exc: ',f16.8,2x,'Evxc: ',f16.8)
 
  9222 format(3x,'Orbital',4x,'Filling',7x,'Eigenvalues',4x,
      $     'Norm(rc->oo)',6x,'Peak')
@@ -903,7 +935,6 @@ subroutine dftsolve(zeff,ixc,exccut,ipsp,ifc,iexit,irel,
          call ghostnl (rvloc)
          call applyaug(rvloc,rvcore(1,iloc))
       endif
-
       call flush(7)
       close(unit=7)
 

diff --git a/src/atom/do_ae.c b/src/atom/do_ae.c
@@ -89,7 +89,10 @@ int do_ae(param_t *param, char *logfile){
     irelxc=0;    
     if (param->ixc < 0) {
       hfsolve_(&param->z,&param->ixc,&exccut_temp,&ipsp,&ifc,&iexit,&irel,&iprint);
-    }else {
+    } else if (param->ixc==7){
+//         printf("z=%d ixc=%d exccut_temp=%lf ipsp=%d ifc=%d iexit=%d irel=%d iprint=%d \n",param->z, param->ixc,exccut_temp,ipsp,iexit, irel,iprint);
+      hfsolve_(&param->z,&param->ixc,&exccut_temp,&ipsp,&ifc,&iexit,&irel,&iprint);
+    } else {
       dftsolve_(&param->z,&param->ixc,&exccut_temp,&ipsp,&ifc,&iexit,&irel,&irelxc,&iprint);
     }
     if (iexit) {
@@ -98,9 +101,13 @@ int do_ae(param_t *param, char *logfile){
     }
 
   } else {
-
+/*   
+ *   Relativity is not implemented for hybrid functionals
+ *   Generating the non-relativistic potentials
+*/      
     fprintf(fp_log," Performing relativistic AE calculation...  \n" );
     fclose(fp_log);
+
 
     /* turn the l-orbitals into j-orbitals */
     relorbae(param,config,logfile);
@@ -117,6 +124,14 @@ int do_ae(param_t *param, char *logfile){
 
     if (param->ixc < 0) {
       dfsolve_(&param->z,&param->ixc,&exccut_temp,&ipsp,&ifc,&iexit,&irel,&iprint);
+    }else if (param->ixc == 7){
+//  non-relativistic calculation for hybrid functionals
+      irel=0;
+      irelxc=0;
+ //     printf("z=%d ixc=%d exccut_temp=%lf ipsp=%d ifc=%d iexit=%d irel=%d iprint=%d \n",param->z, param->ixc,exccut_temp,ipsp,iexit, irel,iprint);
+      nrelorbae(param,config,logfile);
+      startae(param, aorb_.norb);
+      hfsolve_(&param->z,&param->ixc,&exccut_temp,&ipsp,&ifc,&iexit,&irel,&iprint);
     }else {
       dftsolve_(&param->z,&param->ixc,&exccut_temp,&ipsp,&ifc,&iexit,&irel,&irelxc,&iprint);
     }

diff --git a/src/atom/do_nl.c b/src/atom/do_nl.c
@@ -31,7 +31,7 @@
 #include "energy.h"           /* this is for saving the energies */
 
 /* report feature */
-
+//comment by JY hfsolve_ defined here
 static char report[8000];
 void nrelorbnl(param_t *param, int, char *); 
 char * write_reportnl(param_t *param, char *rp,int,double temp_eigen[], double temp_norm[], int, int);
@@ -104,23 +104,26 @@ int do_nl(param_t *param, char *logfile, int insl ){
   irel=0;
   irelxc=0;
   iprint=1;
-  if (param->ixc < 0) {
+  if (param->ixc < 0 || param->ixc == 7) {
     nlpot2_.inl = 0;  
     if (insl != 0 ){
-      printf("!!NOTE!!: All Hartree-Fock tests are done with the semi-local form of the potential \n");
+      printf("!!NOTE!!: All Hartree-Fock and hybrid functional tests are done with the semi-local form of the potential \n");
       printf("!!NOTE!!: Use the 'sl' command instead of the 'nl' command to remove this warning...\n");
 
       fp_log = fopen(logfile, "a");
 
-      fprintf(fp_log,"!!NOTE!!: All Hartree-Fock tests are done with the semi-local form of the potential \n");
+      fprintf(fp_log,"!!NOTE!!: All Hartree-Fock and hybrid functional tests are done with the semi-local form of the potential \n");
       fprintf(fp_log,"!!NOTE!!: Use the 'sl' command instead of the 'nl' command to remove this warning...\n");
 
       fclose(fp_log);
 
     }
     insl=0;
+    //printf("zeff=%lf ixc=%i exccut_temp=%lf ipsp=%i ifc=%i iexit=%i irel=%i iprint=%i",zeff,param->ixc,exccut_temp,ipsp,ifc,iexit,irel,iprint);
     hfsolve_(&zeff,&param->ixc,&exccut_temp,&ipsp,&ifc,&iexit,&irel,&iprint);
   }else {
+    //printf("do nonlocal \n");
+    //printf("zeff=%lf ixc=%i exccut_temp=%lf ipsp=%i ifc=%i iexit=%i irel=%i irelxc=%i iprint=%i",zeff,param->ixc,exccut_temp,ipsp,ifc,iexit,irel,irelxc,iprint);
     dftsolve_(&zeff,&param->ixc,&exccut_temp,&ipsp,&ifc,&iexit,&irel,&irelxc,&iprint);
   }
 
@@ -436,7 +439,7 @@ void writeNL(param_t *param) {
   }
   fclose(fp);
 
-  if (param->ixc >= 0) {
+  if (param->ixc >= 0 && param->ixc != 7) {
     sprintf(filename, "%s.vi_plt", param->name);
     fp = fopen(filename, "a");
     for (j=0;j<param->ngrid;j++){
@@ -509,7 +512,7 @@ void writeSL(param_t *param) {
   }
   fclose(fp);
 
-  if (param->ixc >= 0) {
+  if (param->ixc >= 0 && param->ixc != 7) {
     sprintf(filename, "%s.vi_plt", param->name);
     fp = fopen(filename, "a");
     for (j=0;j<param->ngrid;j++){