Add makefun constructor

devel_tools/makefun_factory/construct.py

@@ -0,0 +1,25 @@
+import glob
+
+# Load header
+with open('makefun_header.f90', 'r') as f:
+    header = f.read()
+
+# Load footer
+with open('makefun_footer.f90', 'r') as f:
+    footer = f.read()
+
+# Load orbitals
+# search for all file that fits name pattern orb_*.f90
+orb_files = glob.glob('orb_*.f90')
+# read files
+orbitals = {f.replace("orb_", "").replace(".f90",""): open(f, 'r').read() for f in orb_files}
+
+# Ensamlble makefun.f90
+with open('makefun_out.f90', 'w') as f:
+    f.write(header)
+    f.write('select case (iopt)\n')
+    for k, v in orbitals.items():
+        f.write(f'case ({k})\n')
+        f.write(v)
+    f.write(footer)
+

devel_tools/makefun_factory/makefun_footer.f90

@@ -0,0 +1,10 @@
+case default
+write(6,*) 'WARNING makefun: orbital',iopt,'not found'
+
+iflagerr=1
+
+end select
+! ** ** ** ** ** ** ** END OF JASTROW ORBITALS ** ** ** ** ** ** ** ** *
+
+return
+end

devel_tools/makefun_factory/makefun_header.f90

@@ -0,0 +1,40 @@
+!TL off
+subroutine makefun(iopt,indt,i0,indtmin,indtm,typec,indpar               &
+           &,indorb,indshell,nelskip,z,dd,zeta,r,rmu,distp               &
+           &,iflagnorm_unused,cr)
+    use constants
+    implicit none
+    integer iopt,indt,i,k,nelskip,indpar,indorbp                         &
+           &,indorb,indshell,indshellp,ic,indtmin,i0                     &
+           &,iflagnorm_unused,indparp,indtm,npower,typec                 &
+           &,ii,jj,kk
+    real*8 z(nelskip,i0:*),dd(*),zeta(*),rmu(3,0:indtm)                  &
+           &,r(0:indtm)                                                  &
+           &,distp(0:indtm,20),peff,fun,fun0,fun2                        &
+           &,rp1,rp2,rp3,rp4,rp5,rp6,rp7,rp8                             &
+           &,dd1,dd2,dd3,dd4,dd5,c,cr,funp,fun2p,funb                    &
+           &,peff2,arg,c0,c1,cost,zv(6),yv(6),xv(6),r2,r4,r6 ! up to i
+
+    integer :: count, multiplicity
+    integer, parameter :: max_power = 20
+    real*8 :: powers(3,0:max_power,0:indtm)
+  !
+  ! indorb are the number of orbitals occupied before calling
+  ! this subroutine
+  !
+  ! indpar is the number of variational parameters used
+  ! before calling this subroutine
+  !
+  ! indshell is the index of the last  occupied orbital
+  ! in the shell, characterized by occupation number iocc(indshell)
+  !
+  ! z(i,indt+4)   contains the laplacian of the orbital i
+  ! z(i,indt+mu) contains the gradient of the orbital i (mu=1,2,3)
+  ! In the following given a radial part of the orbital f(r)
+  ! fun=1/r  d f(r)/d r
+
+  !print *,__FILE__
+  !print *,'makefun: iopt=',iopt
+  !print *,'makefun: i=',i0,' a=',indtmin,' b=',indtm
+  !print *,'makefun: indpar=',indpar,' indorb=',indorb,' indshell=',indshell
+  !print *,'makefun: nelskip=',nelskip

devel_tools/makefun_factory/orb_1.f90

@@ -0,0 +1,36 @@
+  ! s orbital
+  ! 
+  ! - angmom = 0
+  ! - type = Slater
+  ! - normalized = yes
+  ! - angtype = spherical
+  ! - npar = 1
+  ! - multiplicity = 1
+  !
+
+  indshellp=indshell+1
+  dd1=dd(indpar+1)
+  c=dd1*dsqrt(dd1)*0.56418958354775628695d0
+
+  indorbp=indorb+1
+  do k=indtmin,indtm
+    distp(k,1)=c*dexp(-dd1*r(k))
+  end do
+
+  do i=i0,indtm
+    z(indorbp,i)=distp(i,1)
+  end do
+
+  if(typec.ne.1) then
+    fun=-dd1*distp(0,1)
+
+    do i=1,3
+      z(indorbp,indt+i)=fun*rmu(i,0)/r(0)
+    end do
+
+    z(indorbp,indt+4)=(-2.d0*dd1/r(0)+dd1**2)*distp(0,1)
+  end if
+
+  indorb=indorbp
+  indpar=indpar+1
+  indshell=indshellp

devel_tools/makefun_factory/orb_10.f90

@@ -0,0 +1,40 @@
+  ! R(r)=r**2*exp(-z1*r)
+
+
+
+  indshellp=indshell+1
+
+  !        if(iocc(indshellp).eq.1) then
+
+  indorbp=indorb+1
+  dd1=dd(indpar+1)
+  !           if(iflagnorm.gt.2) then
+  !           c=dsqrt((2*dd1)**7/720.d0/pi)/2.d0
+  c=dd1**3.5d0*0.11894160774351807429d0
+  !           endif
+
+  do k=indtmin,indtm
+    distp(k,1)=c*dexp(-dd1*r(k))
+  end do
+
+  do i=i0,indtm
+    z(indorbp,i)=distp(i,1)*r(i)**2
+  end do
+
+  if(typec.ne.1) then
+    fun=(2.d0-dd1*r(0))*distp(0,1)
+    fun2=(2.d0-4*dd1*r(0)+(dd1*r(0))**2)*distp(0,1)
+    do i=1,3
+      z(indorbp,indt+i)=fun*rmu(i,0)
+    end do
+    z(indorbp,indt+4)=2.d0*fun+fun2
+  end if
+
+  indorb=indorbp
+
+  !        endif
+  indpar=indpar+1
+  indshell=indshellp
+
+
+  ! 3s double zeta

devel_tools/makefun_factory/orb_100.f90

@@ -0,0 +1,37 @@
+  !     2s single gaussian
+  !     exp(-dd2*r^2)
+
+
+  dd2=dd(indpar+1)
+
+
+  indorbp=indorb+1
+  indshellp=indshell+1
+  do k=indtmin,indtm
+    distp(k,1)=dexp(-dd2*r(k)*r(k))
+  end do
+
+  !           if(iocc(indshellp).eq.1) then
+  do i=i0,indtm
+    z(indorbp,i)=distp(i,1)
+  end do
+  !           endif
+
+
+  if(typec.ne.1) then
+    fun=-dd2*distp(0,1)*2.d0
+    do i=1,3
+      z(indorbp,indt+i)=fun*rmu(i,0)
+    end do
+
+    z(indorbp,indt+4)=2.d0*dd2*(-3.d0+2.d0*dd2*r(0)**2)*             &
+        distp(0,1)
+    !endif for indt
+  end if
+
+  indpar=indpar+1
+  indshell=indshellp
+  indorb=indorbp
+
+
+

devel_tools/makefun_factory/orb_10000:11000.f90

@@ -0,0 +1,2 @@
+  ! Reserved for dummy orbitals
+

devel_tools/makefun_factory/orb_1000:1099.f90

@@ -0,0 +1,43 @@
+!     s gaussian  r**(2*npower)*exp(-alpha*r**2)
+
+npower=iopt-1000
+
+indorbp=indorb+1
+indshellp=indshell+1
+
+dd2=dd(indpar+1)
+do k=indtmin,indtm
+  distp(k,1)=r(k)**(2*npower)*dexp(-dd2*r(k)**2)
+end do
+
+!        if(iocc(indshellp).eq.1) then
+do i=i0,indtm
+  z(indorbp,i)=distp(i,1)
+end do
+!        endif
+
+
+if(typec.ne.1) then
+
+
+  rp1=r(0)**2
+  fun=(npower-dd2*rp1)*distp(0,1)*2.d0/rp1
+  fun2=(npower*(2.d0*npower-1.d0)-                                   &
+      (1.d0+4.d0*npower)*dd2*rp1+2.d0*(dd2*rp1)**2)*                 &
+      distp(0,1)*2.d0/rp1
+
+  !           if(iocc(indshellp).eq.1) then
+  do i=1,3
+    z(indorbp,indt+i)=rmu(i,0)*fun
+  end do
+  z(indorbp,indt+4)=2.d0*fun+fun2
+  !           endif
+
+
+  !endif for indt
+end if
+
+indpar=indpar+1
+indshell=indshell+1
+indorb=indorbp
+

devel_tools/makefun_factory/orb_101.f90

@@ -0,0 +1,39 @@
+  !     2s without cusp condition
+  !     dd1*( dd3 +exp(-dd2*r^2))
+
+
+  dd2=dd(indpar+1)
+  dd3=dd(indpar+2)
+
+  indorbp=indorb+1
+  indshellp=indshell+1
+  do k=indtmin,indtm
+    distp(k,1)=dexp(-dd2*r(k)*r(k))
+  end do
+
+  !           if(iocc(indshellp).eq.1) then
+  do i=i0,indtm
+    z(indorbp,i)=distp(i,1)+dd3
+  end do
+  !           endif
+
+
+  if(typec.ne.1) then
+    fun=-dd2*distp(0,1)*2.d0
+
+    do i=1,3
+      z(indorbp,indt+i)=fun*rmu(i,0)
+    end do
+
+    z(indorbp,indt+4)=2.d0*dd2*(-3.d0+2.d0*dd2*r(0)**2)*             &
+        distp(0,1)
+
+
+    !endif for indt
+  end if
+
+  indpar=indpar+2
+  indshell=indshellp
+  indorb=indorbp
+
+

devel_tools/makefun_factory/orb_102.f90

@@ -0,0 +1,52 @@
+  !     2s double gaussian with constant
+  !     (dd3+ exp (-dd2 r^2)+dd4*exp(-dd5*r^2))
+
+
+
+  dd2=dd(indpar+1)
+  dd3=dd(indpar+2)
+  dd4=dd(indpar+3)
+  dd5=dd(indpar+4)
+
+  indorbp=indorb+1
+  indshellp=indshell+1
+  do k=indtmin,indtm
+    distp(k,1)=dexp(-dd2*r(k)*r(k))
+    distp(k,2)=dexp(-dd5*r(k)*r(k))
+  end do
+
+  !           if(iocc(indshellp).eq.1) then
+  do i=i0,indtm
+    z(indorbp,i)=(distp(i,1)+dd3+dd4*distp(i,2))
+!              write(6,*) ' function inside = ',z(indorbp,i)
+end do
+!           endif
+
+
+if(typec.ne.1) then
+  fun=-2.d0*(dd2*distp(0,1)+dd5*dd4*distp(0,2))
+  fun2=r(0)**2
+
+  !              write(6,*) ' fun inside = ',fun,fun2
+
+  do i=1,3
+    z(indorbp,indt+i)=fun*rmu(i,0)
+  end do
+
+  z(indorbp,indt+4)=2.d0*(dd2*(-3.d0+2.d0*dd2*fun2)*                 &
+      distp(0,1)+dd5*dd4*(-3.d0+2.d0*dd5*fun2)*distp(0,2))
+
+  !           write(6,*) ' lap 106 =',z(indorbp,indt+4)
+
+
+  !           stop
+
+  !endif for indt
+end if
+
+indpar=indpar+4
+indshell=indshellp
+indorb=indorbp
+
+
+

devel_tools/makefun_factory/orb_103.f90

@@ -0,0 +1,49 @@
+!     2p single gaussian
+
+
+
+dd2=dd(indpar+1)
+
+do k=indtmin,indtm
+  distp(k,1)=dexp(-dd2*r(k)**2)
+end do
+
+!        indorbp=indorb
+
+do ic=1,3
+  !           if(iocc(indshell+ic).eq.1) then
+  indorbp=indorb+ic
+  do i=i0,indtm
+    z(indorbp,i)=rmu(ic,i)*distp(i,1)
+  end do
+  !           endif
+end do
+
+
+if(typec.ne.1) then
+  fun0=distp(0,1)
+  fun=-dd2*distp(0,1)*2.d0
+  fun2=2.d0*dd2*(-1.d0+2.d0*dd2*r(0)**2)*                            &
+      distp(0,1)
+
+  !              indorbp=indorb
+  do ic=1,3
+    !                if(iocc(indshell+ic).eq.1) then
+    indorbp=indorb+ic
+    do i=1,3
+      z(indorbp,indt+i)=rmu(ic,0)*rmu(i,0)*                          &
+          fun
+    end do
+    z(indorbp,indt+ic)=z(indorbp,indt+ic)+fun0
+    z(indorbp,indt+4)=rmu(ic,0)                                      &
+        *(4.d0*fun+fun2)
+    !                 endif
+  end do
+
+  !endif for indt
+end if
+
+indpar=indpar+1
+indshell=indshell+3
+indorb=indorbp
+