lukasplevac/easyconfigs-it4i
1
0
Fork 0
mirror of https://code.it4i.cz/sccs/easyconfigs-it4i.git synced 2025-04-15 19:28:06 +01:00
Code Issues Packages Projects Releases Wiki Activity
easyconfigs-it4i/m/MRCC/mrcc_files/2/combin.f
Lukas Krupcik eb1ea94130 merge master
2024-07-25 10:27:17 +02:00

6201 lines
191 KiB
Fortran
Raw Blame History

************************************************************************
integer*8 function fact(n)
************************************************************************
* This function gives back n factorial *
************************************************************************
implicit none
integer n,i
integer*8 isum
C
isum=1
do i=2,n
isum=isum*i
enddo
fact=isum
return
end
C
************************************************************************
real*8 function factd(n)
************************************************************************
* This function return n factorial in double precision
************************************************************************
implicit none
integer n,i
real*8 isum
C
isum=1.d0
do i=2,n
isum=isum*dble(i)
enddo
factd=isum
return
end
C
************************************************************************
integer*8 function dfact(n)
************************************************************************
* This function gives back n double factorial *
************************************************************************
implicit none
integer n,i
integer*8 isum
C
isum=1
if(mod(n,2).eq.0) then
do i=1,n/2
isum=isum*2*i
enddo
else
do i=1,(n+1)/2
isum=isum*(2*i-1)
enddo
endif
dfact=isum
C
return
end
C
************************************************************************
real*8 function ddfact(n)
************************************************************************
* This function gives back n double factorial in double precision *
************************************************************************
implicit none
integer n,i
real*8 isum
C
isum=1.d0
if(mod(n,2).eq.0) then
do i=1,n/2
isum=isum*2.d0*dble(i)
enddo
else
do i=1,(n+1)/2
isum=isum*(2.d0*dble(i)-1.d0)
enddo
endif
ddfact=isum
C
return
end
C
************************************************************************
integer*8 function dfact2(n)
************************************************************************
* This function gives back (2n-1)!! *
************************************************************************
implicit none
integer n,i
integer*8 isum
C
isum=1
do i=3,2*n-1,2
isum=isum*i
enddo
dfact2=isum
return
end
C
************************************************************************
integer*8 function binom(m,n)
************************************************************************
* This function gives back the binomial coefficient *
************************************************************************
implicit none
integer m,n,i
integer*8 isum,fact
C
isum=1
do i=m-n+1,m
isum=isum*i
enddo
binom=isum/fact(n)
return
end
C
************************************************************************
real*8 function lfact(n)
************************************************************************
* This function gives back the logarithm of n factorial *
************************************************************************
implicit none
integer n,i
real*8 sum
C
sum=0.d0
do i=2,n
sum=sum+dlog(dble(i))
enddo
lfact=sum
return
end
C
************************************************************************
real*8 function lbinom(m,n)
************************************************************************
* This function gives back the logarithm of the binomial coefficient *
************************************************************************
implicit none
integer m,n,i
real*8 sum,lfact
C
sum=0
do i=m-n+1,m
sum=sum+dlog(dble(i))
enddo
lbinom=sum-lfact(n)
return
end
C
************************************************************************
subroutine ifillzero(v,n)
************************************************************************
* This subroutine fills an integer array vectorized with zeros *
************************************************************************
implicit none
integer n,i,v(*)
C
c C$OMP PARALLEL DO
c C$OMP& DEFAULT(SHARED)
do i=1,n
v(i)=0
enddo
c C$OMP END PARALLEL DO
C
return
end
C
************************************************************************
subroutine dfillzero(v,n)
************************************************************************
* This subroutine fills an real*8 array vectorized with zeros *
************************************************************************
implicit none
integer n,i
real*8 v(*)
C
c C$OMP PARALLEL DO
c C$OMP& DEFAULT(SHARED)
do i=1,n
v(i)=0.d0
enddo
c C$OMP END PARALLEL DO
C
return
end
C
************************************************************************
subroutine dfillzer1(v,n)
************************************************************************
* This subroutine fills an real*8 array vectorized with zeros *
************************************************************************
implicit none
integer n,i
real*8 v(*)
C
do i=1,n
v(i)=0.d0
enddo
C
return
end
C
************************************************************************
subroutine ishift(v,n)
************************************************************************
* v=v+1 *
************************************************************************
implicit none
integer n,i,v(*)
C
c C$OMP PARALLEL DO
c C$OMP& DEFAULT(SHARED)
do i=1,n
v(i)=v(i)+1
enddo
c C$OMP END PARALLEL DO
C
return
end
C
************************************************************************
subroutine skip1(ilev,file)
************************************************************************
* This subroutine reads the non-used lines from the formula-file. *
************************************************************************
implicit none
integer n,i,isave,nn,ilev,file
C
ilev=ilev+1
read(file,*) n
do nn=1,n
if(ilev.eq.4) then
read(file,*) i,i,i,i,i,i,i,i,i,i,i,i,isave,i,i,i,i,i,i,i,i
else
read(file,*) i,i,i,i,i,i,i,i,i,i,i,isave,i,i,i,i,i,i,i,i
if(isave.gt.0) call skip2(ilev,file)
endif
enddo
ilev=ilev-1
C
return
end
C
************************************************************************
subroutine skip2(ilev,file)
************************************************************************
* This subroutine reads the non-used lines from the formula-file. *
************************************************************************
implicit none
integer ilev,file
C
call skip1(ilev,file)
C
return
end
C
************************************************************************
subroutine comb(n,nact,icomb,ilev,in,inact,ip,ncomb,ic)
************************************************************************
* Generate combinations *
************************************************************************
implicit none
integer n,nact,icomb(0:n,*),ilev,in,inact,ncomb,ic(*),ip
C
call ifillzero(ic,n)
ilev=0
ncomb=0
call comb1(n,nact,icomb,ilev,in,inact,ip,ncomb,ic)
C
return
end
C
************************************************************************
subroutine comb1(n,nact,icomb,ilev,in,inact,ip,ncomb,ic)
************************************************************************
* Generate combinations *
************************************************************************
implicit none
integer n,nact,icomb(0:n,*),ilev,in,inact,ncomb,ic(*),i,j,ip
integer is(n),ina
logical log
C
if(ilev.eq.in) then
ina=0
do i=1,nact ! Only for virtual strings (actives first)!!!
if(ic(i).eq.1) ina=ina+1
enddo
if(ina.ne.inact) return
C
ncomb=ncomb+1
do i=1,n
icomb(i,ncomb)=ic(i)
is(i)=ic(i)
enddo
icomb(0,ncomb)=1
log=.true.
do while(log)
log=.false.
do i=1,n-1
if(is(i).lt.is(i+1)) then
icomb(0,ncomb)=-icomb(0,ncomb)
log=.true.
j=is(i)
is(i)=is(i+1)
is(i+1)=j
endif
enddo
enddo
return
endif
C
ilev=ilev+1
do i=ip+1,n
ic(i)=1
call comb2(n,nact,icomb,ilev,in,inact,i,ncomb,ic)
ic(i)=0
enddo
ilev=ilev-1
C
return
end
C
************************************************************************
subroutine comb2(n,nact,icomb,ilev,in,inact,ip,ncomb,ic)
************************************************************************
* Generate combinations *
************************************************************************
implicit none
integer n,nact,icomb(0:n,*),ilev,in,inact,ncomb,ic(*),ip
C
call comb1(n,nact,icomb,ilev,in,inact,ip,ncomb,ic)
C
return
end
C
#if defined (OMP)
************************************************************************
subroutine ompred(w,n)
************************************************************************
* "Manual" reduction for OpenMP
************************************************************************
#include "MRCCCOMMON"
integer n,i
real*8 w(n,xyzsiz1)
C
do i=2,xyzsiz1
call daxpy(n,1.d0,w(1,i),1,w(1,1),1)
enddo
C
return
end
C
#endif
************************************************************************
subroutine trimloc(beg,end)
************************************************************************
* Local trim function
************************************************************************
implicit none
integer lb,le,i
character*1 beg(*),end(*)
C
lb=0
le=0
do i=1,128
if(beg(i).ne.' ') lb=i
if(end(i).ne.' ') le=i
enddo
do i=lb+1,min(128,lb+le)
beg(i)=end(i-lb)
enddo
do i=lb+le+1,128
beg(i)=' '
enddo
C
return
end
C
************************************************************************
integer function intadd(n)
************************************************************************
* This function returns the address of dcore(n) wrt icore *
************************************************************************
implicit none
integer n,iintfp
#if defined (INT64)
parameter(iintfp=1)
#else
parameter(iintfp=2)
#endif
C
intadd=iintfp*(n-1)+1
C
return
end
C
************************************************************************
integer function dbladd(n)
************************************************************************
* This function returns the address of icore(n) wrt dcore *
************************************************************************
#include "MRCCCOMMON"
integer n,i
C
i=mod(n-1,iintfp)
if(i.ne.0) i=iintfp-i
dbladd=(n-1+i)/iintfp+1
C
return
end
C
************************************************************************
subroutine putint(ifile,irec,vec,ilength)
************************************************************************
* Write an integer array to a direct access file *
* ifile - unit number *
* irec - record number in direct access file *
* vec - array to be read *
* ilength - length of the array *
************************************************************************
#include "MRCCCOMMON"
integer ifile,irec,ilength,nrec,i,j,jrec,vec(*)
C
if(ilength.eq.0) return
j=mod(ilength,iibufln)
nrec=(ilength-j)/iibufln-1
do jrec=0,nrec
write(ifile,rec=irec+jrec) (vec(jrec*iibufln+i),i=1,iibufln)
enddo
if(j.ne.0) write(ifile,rec=irec+nrec+1)
$(vec((nrec+1)*iibufln+i),i=1,j)
C
return
end
C
************************************************************************
subroutine getint(ifile,irec,vec,ilength)
************************************************************************
* Read an integer array from a direct access file *
* ifile - unit number *
* irec - record number in direct access file *
* vec - array to be read *
* ilength - length of the array *
************************************************************************
#include "MRCCCOMMON"
integer ifile,irec,ilength,nrec,i,j,jrec,vec(*)
C
if(ilength.eq.0) return
j=mod(ilength,iibufln)
nrec=(ilength-j)/iibufln-1
do jrec=0,nrec
read(ifile,rec=irec+jrec) (vec(jrec*iibufln+i),i=1,iibufln)
enddo
if(j.ne.0) read(ifile,rec=irec+nrec+1)
$(vec((nrec+1)*iibufln+i),i=1,j)
C
return
end
C
************************************************************************
subroutine putlst(ifile,irec,vec,ilength)
************************************************************************
* Write a double precision array to a direct access file *
* ifile - unit number *
* irec - record number in direct access file *
* vec - array to be read *
* ilength - length of the array *
************************************************************************
#include "MRCCCOMMON"
integer ifile,irec,ilength,nrec,i,j,jrec
real*8 vec(*)
C
if(ilength.eq.0) return
j=mod(ilength,ibufln)
nrec=(ilength-j)/ibufln-1
do jrec=0,nrec
write(ifile,rec=irec+jrec) (vec(jrec*ibufln+i),i=1,ibufln)
enddo
if(j.ne.0) write(ifile,rec=irec+nrec+1)
$(vec((nrec+1)*ibufln+i),i=1,j)
C
return
end
C
************************************************************************
subroutine getlst(ifile,irec,vec,ilength)
************************************************************************
* Read a double precision array from a direct access file *
* ifile - unit number *
* irec - record number in direct access file *
* vec - array to be read *
* ilength - length of the array *
************************************************************************
#include "MRCCCOMMON"
integer ifile,irec,ilength,nrec,i,j,jrec
real*8 vec(*)
C
if(ilength.eq.0) return
j=mod(ilength,ibufln)
nrec=(ilength-j)/ibufln-1
do jrec=0,nrec
read(ifile,rec=irec+jrec) (vec(jrec*ibufln+i),i=1,ibufln)
enddo
if(j.ne.0) read(ifile,rec=irec+nrec+1)
$(vec((nrec+1)*ibufln+i),i=1,j)
C
return
end
C
************************************************************************
subroutine getput(ifile,irec,vec,ilength)
************************************************************************
* Adds vec and a vector from ifile together and writes back into ifile
************************************************************************
#include "MRCCCOMMON"
integer ifile,irec,ilength,nrec,i,j
real*8 vec(*)
C
if(ilength.eq.0) return
nrec=0
j=irec
do while (nrec.lt.ilength)
read(ifile,rec=j) ibuf
c C$OMP PARALLEL DO
c C$OMP& DEFAULT(SHARED)
do i=1,min(ibufln,ilength-nrec)
ibuf(i)=ibuf(i)+vec(nrec+i)
enddo
c C$OMP END PARALLEL DO
write(ifile,rec=j) ibuf
nrec=nrec+ibufln
j=j+1
enddo
C
return
end
C
************************************************************************
subroutine getkeym(key,nlk,val,nlv)
************************************************************************
* Get the value of a keyword from inputfile
************************************************************************
#include "MRCCCOMMON"
character*1 key(nlk),val(nlv),line(512)
character*4 line4
character*8 line8
character*512 line1
integer nlk,nlv,i,n
logical ll
equivalence(line,line1,line4,line8)
C
val(1:nlv)=' '
line(1:512)=' '
rewind(minpfile)
do
read(minpfile,'(512a1)',end=1000) line
call lowercase(line,line,512)
line1=adjustl(line1)
if(line4.eq.'geom') then
line(5)='='
else if(line8.eq.'pointcha') then
return
endif
n=1
do while(line(n).ne.'='.and.line(n).ne.' '.and.n.lt.512)
n=n+1
enddo
if(n-1.eq.nlk) then
ll=.true.
do i=1,nlk
ll=ll.and.line(i).eq.key(i)
enddo
if(ll) then
i=1
do while(line(i).ne.'='.and.i.lt.512)
line(i)=' '
i=i+1
enddo
line(i)=' '
line1=adjustl(line1)
i=1
do while(line(i).ne.' ')
val(i)=line(i)
i=i+1
enddo
return
endif
endif
enddo
1000 continue
C
return
end
C
************************************************************************
subroutine getkey(key,nlk,val,nlv)
************************************************************************
* Get the value of a keyword from keyword file
************************************************************************
#include "MRCCCOMMON"
character*1 key(nlk),val(nlv),line(512)
character*512 line1
integer nlk,nlv,i,n
logical ll
equivalence(line,line1)
C
open(keywdfile,file='KEYWD',status='old')
rewind(keywdfile)
val(1:nlv)=' '
line(1:512)=' '
do
read(keywdfile,'(512a1)',end=1000) line
call lowercase(line,line,512)
line1=adjustl(line1)
n=1
do while(line(n).ne.'='.and.line(n).ne.' '.and.n.le.nlk)
n=n+1
enddo
if(n-1.eq.nlk) then
ll=.true.
do i=1,nlk
ll=ll.and.line(i).eq.key(i)
enddo
if(ll) then
i=1
do while(line(i).ne.'='.and.i.lt.512)
line(i)=' '
i=i+1
enddo
line(i)=' '
line1=adjustl(line1)
i=1
do while(line(i).ne.' '.and.i.le.nlv)
val(i)=line(i)
i=i+1
enddo
goto 1000
endif
endif
enddo
1000 continue
close(keywdfile)
C
return
end
************************************************************************
double precision function dgetkey(key,nlk,nlv) result(dkey)
************************************************************************
* Reads a keyword and converts to double precision
************************************************************************
implicit none
integer nlk,nlv
character*1 key(nlk)
character*128 c
c=' '
call getkey(key,nlk,c,nlv)
read(c,*) dkey
end function
C
************************************************************************
subroutine unknown(key,nlen)
************************************************************************
* Exit if the value of a keyword is not correct
************************************************************************
#include "MRCCCOMMON"
integer nlen,i
character*1 key(nlen)
C
write(iout,*)
write(iout,*)
$'Unknown option for keyword "',(key(i),i=1,nlen),'"!'
call mrccend(1)
C
return
end
C
************************************************************************
subroutine changekey(key,nlk,val,nlv)
************************************************************************
* Change a keyword in keyword file
************************************************************************
#include "MRCCCOMMON"
character*1 key(nlk),val(nlv),line(512)
character*512 line1
integer nlk,nlv,i,n
logical ll
equivalence(line,line1)
C
open(keywdfile,file='KEYWD',status='old')
open(scrfile7,status='unknown')
rewind(keywdfile)
rewind(scrfile7)
do i=1,nlk
line(i)=key(i)
enddo
do
read(keywdfile,'(512a1)',end=1000) line
line1=adjustl(line1)
n=1
do while(line(n).ne.'='.and.line(n).ne.' ')
n=n+1
enddo
if(n-1.eq.nlk) then
ll=.true.
do i=1,nlk
ll=ll.and.line(i).eq.key(i)
enddo
if(ll) then
do i=1,nlv
line(n+i)=val(i)
enddo
do i=n+nlv+1,512
line(i)=' '
enddo
endif
endif
write(scrfile7,'(512a1)') line
enddo
1000 continue
rewind(keywdfile)
rewind(scrfile7)
do
read(scrfile7,'(512a1)',end=1001) line
line1=adjustl(line1)
write(keywdfile,*) trim(line1)
enddo
1001 continue
close(keywdfile)
close(scrfile7,status='delete')
C
return
end
C
************************************************************************
subroutine lowercase(ch,chl,nlen)
************************************************************************
* Convert upper case to lower case
************************************************************************
implicit none
integer nlen,i,ii
character*1 ch(*),chl(*)
C
do i=1,nlen
chl(i)=ch(i)
ii=ichar(chl(i))
if(ii.ge.65.and.ii.le.90) chl(i)=char(ii+32)
enddo
C
return
end
C
************************************************************************
subroutine uppercase(ch,chl,nlen)
************************************************************************
* Convert upper case to lower case
************************************************************************
implicit none
integer nlen,i,ii
character*1 ch(*),chl(*)
C
do i=1,nlen
chl(i)=ch(i)
ii=ichar(chl(i))
if(ii.ge.97.and.ii.le.122) chl(i)=char(ii-32)
enddo
C
return
end
C
#if defined (MPI)
************************************************************************
subroutine mrccini_mpi
************************************************************************
* Initiates MPI
************************************************************************
#include "MRCCCOMMON"
character*16 devnul
common/output/ devnul
integer bcast_size, provided
integer :: mpi_err
integer mpi_rank, mpi_size, nlen, parent
logical mpi_in
character*8 rank_ch
character*512 wdir_parent
character*(MPI_MAX_PROCESSOR_NAME) pname
character*21 prov_ch
logical dodiis
common /pd/ dodiis
tim=0.d0
wct=0.d0
call mtime(tim,wct)
timold=tim
wctold=wct
timvold=tim
wctvold=wct
! iout = 6
mem_allocated = .false.
! del_wdir = .false.
! master_thread = .true.
! master_rank = 0
! inbcast = .true.
bcast_comm = MPI_COMM_NULL
call signalinitchild ! initialise response to SIG's
call MPI_Initialized(mpi_in, mpi_err)
if(.not. mpi_in) then
call MPI_Init_thread(MPI_THREAD_SERIALIZED, !MULTIPLE,
$ provided, mpi_err)
if(provided .ne. MPI_THREAD_SERIALIZED) then !MULTIPLE) then
if(provided .eq. MPI_THREAD_SINGLE) then
prov_ch = 'MPI_THREAD_SINGLE '
elseif(provided .eq. MPI_THREAD_FUNNELED) then
prov_ch = 'MPI_THREAD_FUNNELED '
elseif(provided .eq. MPI_THREAD_SERIALIZED) then
prov_ch = 'MPI_THREAD_SERIALIZED'
elseif(provided .eq. MPI_THREAD_MULTIPLE) then
prov_ch = 'MPI_THREAD_MULTIPLE '
else
prov_ch = 'unknown '
end if
write(6, '("WARNING: MPI initialized at ",
$ a, " thread support level.")') trim(prov_ch)
end if
end if
call MPI_Comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_err)
call MPI_Comm_size(MPI_COMM_WORLD, mpi_size, mpi_err)
call MPI_Get_processor_name(pname, nlen, mpi_err)
call MPI_Comm_get_parent(parent, mpi_err)
call getcwd(wdir)
if(parent .eq. MPI_COMM_NULL) then
if(mpi_rank .eq. 0) then
output_folder = 'mpi_output'
wdir_parent = wdir
call system('mkdir -p ' // trim(output_folder) // char(0))
end if
parent = MPI_COMM_WORLD
end if
c if(parent .ne. MPI_COMM_NULL) then
call MPI_Bcast(wdir_parent, 512, MPI_CHARACTER, 0, parent,
$ mpi_err)
if(trim(wdir) .ne. trim(wdir_parent)) then
call system('mkdir -p ' // trim(wdir_parent) //
$ ' > /dev/null')
call chdir(trim(wdir_parent))
call getcwd(wdir)
if(wdir_parent .ne. wdir) then
write(*,
$ '(" Error: Couldn''t create working directory on ",a)'),
$ trim(pname)
call mrccend(1)
end if
c end if
c if(.not. master_thread) del_wdir = .true.
end if
call MPI_Bcast(output_folder, 32, MPI_CHARACTER, 0, parent,
$ mpi_err)
if(mpi_size .gt. 1) then
call choose_master_process
call create_bcast_group(MPI_COMM_WORLD,master_rank,bcast_comm,
$ inbcast,0)
del_wdir = inbcast .and. .not.master_thread
else
master_rank = 0
inbcast = .true.
master_thread = .true.
del_wdir = .false.
end if
dodiis = .false. !inbcast
c if(inbcast .and. mpi_size .gt. 1) then
c call MPI_Comm_size(bcast_comm, bcast_size, mpi_err)
c if(bcast_size .ne. mpi_size) dodiis = .false.
c end if
! redirect output
write(rank_ch, '(i8)') mpi_rank
if(.not. master_thread) then
outfilename = 'mrcc.out.'//trim(adjustl(rank_ch))
iout = outfile
open(iout, file = outfilename, access = 'append')
devnul = '>>' // trim(outfilename)
else
outfilename = ' '
devnul = ' '
iout = 6
tag = -1
end if
tag = 0
if(mpi_size .ne. 1) then
if(master_thread)
$ write(iout, '(" MPI parallel version",/)')
call MPI_Barrier(MPI_COMM_WORLD, mpi_err)
call sleep_f(1)
write(rank_ch, '(i8)') mpi_rank
if(master_thread) then
write(*,"(' rank ', a4, ' running on ', a20,
$ ' master')")
$ trim(adjustl(rank_ch)), pname
else
write(*,"(' rank ', a4, ' running on ', a20)")
$ trim(adjustl(rank_ch)), pname
end if
call MPI_Barrier(MPI_COMM_WORLD, mpi_err)
call sleep_f(1)
if(master_thread) write(iout, *)
end if
if(master_thread) then
open(unit=exitfile,file='EXIT',status='unknown')
write(exitfile,"(i5,2x,2(1pe15.8))") 1,0.d0,0.d0
close(exitfile)
end if
call addPidToList
end subroutine
************************************************************************
subroutine choose_master_process
************************************************************************
#include "MRCCCOMMON"
integer mpi_rank, mpi_size, mpi_err, nlen, parent, i
character*(MPI_MAX_PROCESSOR_NAME) pname
logical mpi_in, file_exists
logical, allocatable :: master_node(:)
character(len=MPI_MAX_PROCESSOR_NAME) :: id_ref
character*8 date
character*10 time
character*5 zone
integer dt(8)
call MPI_Comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_err)
call MPI_Comm_size(MPI_COMM_WORLD, mpi_size, mpi_err)
call MPI_Get_processor_name(pname, nlen, mpi_err)
call MPI_Comm_get_parent(parent, mpi_err)
allocate(master_node(0:mpi_size-1))
master_node = .false.
if(parent .ne. MPI_COMM_NULL) then
! receive id from dmrcc
call MPI_Bcast(id_ref, MPI_MAX_PROCESSOR_NAME, MPI_CHARACTER,
$ 0, parent, mpi_err)
! choose master process
if(trim(adjustl(pname)) .eq. trim(adjustl(id_ref)))
$ master_node(mpi_rank) = .true.
else
if(mpi_size .gt. 1 .and. mpi_rank .eq. 0) then
write(*, '(/," WARNING: no parent process")')
write(*, '( " assigning rank 0 as master",/)')
end if
if(mpi_rank .eq. 0) then
master_node(0) = .true.
call date_and_time(date, time, zone, dt)
id_ref = pname
end if
call MPI_Bcast(id_ref, MPI_MAX_PROCESSOR_NAME, MPI_CHARACTER,
$ 0, MPI_COMM_WORLD, mpi_err)
end if
call MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL,
$ master_node, 1, MPI_LOGICAL, MPI_COMM_WORLD, mpi_err)
do i = 0, mpi_size-1
if(master_node(i)) exit
end do
master_rank = i
if(i .eq. mpi_size) then
write(*,*) "ERROR: input files not found"
call mrccend(1)
else if(i .eq. mpi_rank) then
master_thread = .true.
else
master_thread = .false.
end if
deallocate(master_node)
end subroutine
************************************************************************
subroutine create_bcast_group(communicator, master_rank,
$ bcast_comm, inbcast, group_id)
************************************************************************
implicit none
#include "mpif.h"
integer mpi_group_world, bcast_group, group_id
integer group_size, communicator, master_rank, bcast_comm
integer mpi_rank, mpi_size, mpi_err, nlen, parent, i
integer, allocatable :: group_members(:)
logical mpi_in, file_exists, ingroup, inbcast
character*8 rank_ch,group_id_c
character*64 id, id_ref, group_idfile
character*(MPI_MAX_PROCESSOR_NAME) pname
data id /''/
character*8 date
character*10 time
character*5 zone
integer dt(8)
integer, parameter :: idfile = 961
! choose one process for each disk
call MPI_Comm_rank(communicator, mpi_rank, mpi_err)
call MPI_Comm_size(communicator, mpi_size, mpi_err)
call MPI_Get_processor_name(pname, nlen, mpi_err)
allocate(group_members(mpi_size))
write(group_id_c,'(i6)') group_id
group_idfile=trim('idfile') // '.' // trim(group_id_c)
c write(*,*) 'rank,idf',mpi_rank, trim(group_idfile)
if(mpi_rank .eq. master_rank) then
call date_and_time(date, time, zone, dt)
id_ref = date // time // zone // pname
open(idfile, file = trim(group_idfile))
write(idfile, *) trim(id_ref)
close(idfile)
inbcast = .true.
end if
call MPI_Bcast(id_ref, 64, MPI_CHARACTER, master_rank,
$ communicator, mpi_err)
group_size = 1
group_members(group_size) = master_rank
call MPI_Barrier(communicator, mpi_err)
do i = 0, mpi_size-1
if(i .eq. master_rank) cycle
if(i .eq. mpi_rank) then
id = ''
inquire(file = trim(group_idfile) , exist = file_exists)
if(file_exists) then
open(idfile, file = trim(group_idfile) )
read(idfile, *) id
close(idfile)
end if
if(id .ne. id_ref) then
open(idfile, file = trim(group_idfile) )
write(idfile, *) trim(id_ref)
close(idfile)
ingroup = .true.
else
ingroup = .false.
end if
inbcast = ingroup
end if
call MPI_Bcast(ingroup, 1, MPI_LOGICAL, i,
$ communicator, mpi_err)
if(ingroup) then
group_size = group_size + 1
group_members(group_size) = i
end if
end do
call MPI_Comm_group(communicator, mpi_group_world, mpi_err)
call MPI_Group_incl(mpi_group_world, group_size, group_members,
$ bcast_group, mpi_err)
call MPI_Comm_create(communicator, bcast_group, bcast_comm,
$ mpi_err)
deallocate(group_members)
c call delete_file('dmrcc.id', inbcast)
call delete_file(trim(group_idfile), inbcast)
end subroutine
************************************************************************
subroutine bcast_file(filename, comm)
************************************************************************
#include "MRCCCOMMON"
integer :: comm
character(len = *) :: filename
integer :: mpi_rank, mpi_err, comm_size
integer :: msg_size
integer :: iostatus
integer, parameter :: buff_size = 1048576 ! 1 MB
character(len = :), allocatable :: buffer
logical :: flag, flag2
logical :: empty, fexists
logical fexist
c integer fsize
if(.not.inbcast) return
call MPI_Initialized(flag, mpi_err)
call MPI_Finalized(flag2, mpi_err)
if(.not. flag .or. flag2) return
if(comm .eq. MPI_COMM_NULL) return
c fpos = 1
call MPI_Comm_size(comm, comm_size, mpi_err)
if(comm_size .eq. 1) return
call MPI_Comm_rank(comm, mpi_rank, mpi_err)
inquire(file = filename, exist = fexist)
allocate(character(len = buff_size) :: buffer)
if(mpi_rank .eq. 0) then
c msg_size = min(buff_size, fsize)
msg_size = buff_size
inquire(file = filename, exist = fexists)
call MPI_Bcast(fexists, 1, MPI_LOGICAL,
$ 0, comm, mpi_err)
open(unit = cpfile, file = filename, access = 'stream')
do
! read and bcast filename
buffer = repeat(' ', buff_size)
read(cpfile, iostat = iostatus) buffer
c read(cpfile, pos = fpos, iostat = iostatus) buffer
c fpos = fpos + buff_size
! end of file
if(iostatus .ne. 0) msg_size = len_trim(buffer)
if(msg_size .eq. 0) exit ! empty or end of file
call MPI_Bcast(msg_size, 1, MPI_INTEGER_MRCC,
$ 0, comm, mpi_err)
call MPI_Bcast(buffer, msg_size, MPI_CHARACTER,
$ 0, comm, mpi_err)
if(iostatus .ne. 0) exit
end do
call MPI_Bcast(0, 1, MPI_INTEGER_MRCC,
$ 0, comm, mpi_err)
c empty = .not.fexists
else
call delete_file(filename, .true.)
close(cpfile)
open(unit = cpfile, file = filename, access = 'stream')
call MPI_Bcast(fexists, 1, MPI_LOGICAL,
$ 0, comm, mpi_err)
c empty = .true.
do
! receive bcast
call MPI_Bcast(msg_size, 1, MPI_INTEGER_MRCC,
$ 0, comm, mpi_err)
! end of file
if(msg_size .eq. 0) exit
c empty = .false.
call MPI_Bcast(buffer, msg_size, MPI_CHARACTER,
$ 0, comm, mpi_err)
! write to file
write(cpfile) buffer(1:msg_size)
c write(cpfile, pos = fpos) buffer(1:msg_size)
c fpos = fpos + msg_size
end do
end if
empty = .not.fexists
close(cpfile)
call delete_file(filename, empty)
deallocate(buffer)
call MPI_Barrier(comm, mpi_err)
end subroutine
************************************************************************
subroutine gather_file(filename, comm)
************************************************************************
c Gathers output files from other ranks. Only works with comm=bcast_comm
#include "MRCCCOMMON"
integer comm
character(len = *) filename
integer :: mpi_err, comm_size, msg_tag, mpi_rank
integer :: msg_size
integer :: iostatus, i
integer :: stat(MPI_STATUS_SIZE)
integer, parameter :: buff_size = 1048576 ! 1 MB
character(len = 4) i_ch
character(len = :), allocatable :: buffer
logical :: flag, flag2
call MPI_Initialized(flag, mpi_err)
call MPI_Finalized(flag2, mpi_err)
if(.not. flag .or. flag2) return
if(comm .eq. MPI_COMM_NULL) return
call MPI_Comm_size(comm, comm_size, mpi_err)
if(comm_size .eq. 1) return
call MPI_Comm_rank(comm, mpi_rank, mpi_err)
allocate(character(len = buff_size) :: buffer)
c if(.not. master_thread) then
if(mpi_rank .ne. master_rank) then
open(unit = cpfile, file = filename, access = 'stream')
msg_tag = 0
do
! read file and send it to master
buffer = repeat(' ', buff_size)
read(cpfile, iostat = iostatus) buffer
if(iostatus .ne. 0) then
! end of file
msg_size = len_trim(buffer)
msg_tag = 1
end if
call MPI_Send(msg_size, 1, MPI_INTEGER_MRCC,
$ master_rank, msg_tag, comm, mpi_err)
call MPI_Send(buffer, msg_size, MPI_CHARACTER,
$ master_rank, msg_tag, comm, mpi_err)
if(msg_tag .ne. 0) exit
end do
close(cpfile)
else
do i = 0, comm_size-1
if(i .eq. master_rank) cycle
write(i_ch, '(i4)') i
open(unit = cpfile, file = trim(adjustl(output_folder)) //
$ '/' // trim(adjustl(filename)) // '.' //
$ trim(adjustl(i_ch)),
$ access = 'stream', position = 'append')
do
! receive bcast
call MPI_Recv(msg_size, 1, MPI_INTEGER_MRCC,
$ i, MPI_ANY_TAG, comm, stat, mpi_err)
call MPI_Recv(buffer, msg_size, MPI_CHARACTER,
$ i, MPI_ANY_TAG, comm, stat, mpi_err)
! write to file
write(cpfile) buffer(1:msg_size)
! end of file
if(stat(MPI_TAG) .ne. 0) exit
end do
close(cpfile)
end do
end if
deallocate(buffer)
end subroutine
************************************************************************
subroutine mpifinalize(ierr)
************************************************************************
* Finalizes MPI
************************************************************************
#include "MRCCCOMMON"
integer ierr
integer mpi_err, request, parent
integer system_error, nlen
integer system
logical mpi_in, mpi_fin
character(len = MPI_MAX_PROCESSOR_NAME) pname
call MPI_Initialized(mpi_in, mpi_err)
call MPI_Finalized(mpi_fin, mpi_err)
if(mpi_in .and. .not.mpi_fin) then
c if(tag.ge.0) then
call MPI_Get_processor_name(pname, nlen, mpi_err)
if(del_wdir) then
call delete_file('pids', inbcast)
call chdir('..')
system_error = system('rmdir ' // trim(wdir))! // ' 2>/dev/null') !!! wdir should have been set in mpiinit
if(system_error .ne. 0) then
write(*,
$ '(" Error: Couldn''t delete working directory on ", a)')
$ trim(pname)
end if
end if
flush(iout)
if(.not. master_thread) close(iout)
if(ierr .eq. -1) call exit(abs(ierr)) ! parent signalled child, no need to notify parent
call MPI_Comm_get_parent(parent, mpi_err)
! signal parents
if(parent .ne. MPI_COMM_NULL) then
call MPI_Ibarrier(parent, request, mpi_err)
call MPI_Wait(request, MPI_STATUS_IGNORE, mpi_err)
call sleep_f(1)
call MPI_Comm_disconnect(parent, mpi_err)
c call MPI_Comm_free(parent, mpi_err)
end if
call MPI_Finalize(mpi_err)
end if
end subroutine
#endif
************************************************************************
subroutine delete_file(filename, delete)
************************************************************************
#include "MRCCCOMMON"
character(len = *) :: filename
logical :: delete
#if defined(MPI)
if(.not.delete) return
#endif
c$OMP CRITICAL
open(delfile, file = trim(filename))
close(delfile, status = 'delete')
c$OMP END CRITICAL
end subroutine
************************************************************************
subroutine file_copy_fortran(origfilename,copiedfilename)
************************************************************************
#include "MRCCCOMMON"
character(len = *) :: origfilename,copiedfilename
c character*4 origfilename,copiedfilename
integer, parameter :: buff_size = 1048576 ! 1 MB
integer :: iostatus,write_size
character(len = :), allocatable :: buffer
logical lll
call delete_file(trim(copiedfilename), .true.)
open(unit = delfile, file = trim(origfilename),
$access = 'stream',IOSTAT=iostatus)
open(unit=cpfile,file=trim(copiedfilename),access='stream')
write_size=buff_size
allocate(character(len = buff_size) :: buffer)
do
buffer = repeat(' ', buff_size)
read(delfile, iostat = iostatus) buffer
if(iostatus.ne.0) write_size = len_trim(buffer)
if(write_size.eq.0) exit ! empty file
write(cpfile) buffer(1:write_size)
if(iostatus.ne.0) exit ! end of file
end do
close(delfile)
close(cpfile)
deallocate(buffer)
end subroutine
************************************************************************
subroutine mrccini
************************************************************************
* Initiates mrcc
************************************************************************
#include "MRCCCOMMON"
C
tim=0.d0
wct=0.d0
call mtime(tim,wct)
timold=tim
wctold=wct
timvold=tim
wctvold=wct
call addPidToList
iout = 6
mem_allocated = .false.
del_wdir = .false.
master_thread = .true.
master_rank = 0
inbcast = .true.
bcast_comm = -1
open(unit=exitfile,file='EXIT',status='unknown')
write(exitfile,"(i5,2x,2(1pe15.8))") 1,0.d0,0.d0
close(exitfile)
call signalinitchild ! initialise response to SIG's
C
return
end
C
************************************************************************
subroutine mrccini_from_subroutine
************************************************************************
* Initialize if mrccini is not executed when calling a program as a subroutine
************************************************************************
#include "MRCCCOMMON"
C
call mtime(tim,wct)
timold=tim
wctold=wct
C
return
end
c
************************************************************************
subroutine mrccend(ierr)
************************************************************************
* Terminates mrcc
************************************************************************
#include "MRCCCOMMON"
integer ierr
C
tim=1.d0
wct=1.d0
call mtime(tim,wct)
call removePidFromList
#if defined(MPI)
if(master_thread) then
#endif
open(unit=exitfile,file='EXIT',status='unknown')
write(exitfile,"(i5,2x,2(1pe15.8))") ierr,tim-timvold,wct-wctvold
close(exitfile)
#if defined(MPI)
end if
call mpifinalize(ierr)
#endif
if(ierr .eq. -2) return ! return to signal handler to call previous handler
call exit(ierr)
C
return
end
C
************************************************************************
subroutine addPidToList
************************************************************************
* Adds PID of the calling process to the list of running processes
* stored in file pids
************************************************************************
#include "MRCCCOMMON"
integer :: getpid_f
open(unit = pidfile, file = 'pids', access = 'append')
write(pidfile, *) getpid_f()
close(unit = pidfile)
end subroutine
C
************************************************************************
subroutine removePidFromList
************************************************************************
* Removes PID of the calling process from the list of running processes
* stored in file pids
************************************************************************
#include "MRCCCOMMON"
integer :: ownPid, getpid_f
character(len = 8) :: pidChar
ownPid = getpid_f()
write(pidChar, '(i8)') ownPid
call system('sed -i "/' // trim(adjustl(pidChar)) //
$ '/d" pids 2> /dev/null')
end subroutine
C
************************************************************************
subroutine timer
************************************************************************
* Returns CPU time minutes *
************************************************************************
#include "MRCCCOMMON"
C
tim=1.d0
wct=1.d0
call mtime(tim,wct)
write(iout,"(' CPU time [min]:',f10.3,
$' Wall time [min]:',f10.3)") tim/60.d0,wct/60.d0
C
return
end
C
************************************************************************
subroutine mtime(tt,ww)
************************************************************************
* Machine-dependent TIME routine *
* Returns CPU and wall time in seconds via double precision tt and wct *
************************************************************************
implicit none
integer wzero,w,ps
real*8 tt,tzero,t,ww
save tzero
save wzero
c real*4 t
C# AIX
C real*8 timef
C# LINUX
c call CPU_TIME(t)
c tt=t
call CPU_TIME(t)
if(tt.eq.0.d0) tzero=t
tt=t-tzero
call SYSTEM_CLOCK(w,ps)
if(ww.eq.0.d0) wzero=w
ww=dble(w-wzero)/dble(ps)
C# AIX
C tt=timef()/1000.d0
return
end
C
************************************************************************
subroutine ishell(string)
************************************************************************
* Call external programs *
************************************************************************
#include "MRCCCOMMON"
character*(*) string
integer*4 istat
#if !defined (gfortran) && !defined (G95)
integer*4 system
external system
#endif
#if defined (Solaris10) || defined (gfortran)
character cpath*4096, cmd*5000
if( string(1:5).eq.'echo '
$ .or. string(1:4).eq.'sed '
$ .or. string(1:4).eq.'cat ') then
cmd=string
goto 100
endif
call getenv('PATH', cpath)
C Test if cpath was long enough to contain PATH
if(len_trim(cpath).gt.len(cpath)-5) then
write(iout,*)
write(iout,*) 'ishell: cpath variable too short.'
write(iout,*) 'Program will stop.'
call flush(iout)
call exit(1)
endif
C Build the tricky command
cmd="bash -c "//"'"//"export PATH="//cpath(1:len_trim(cpath))
$ //"; "//string(1:len_trim(string))//"'"
C Test if cmd was long enough to expand PATH and string
if(len_trim(cmd).gt.len(cmd)-5) then
write(iout,*)
write(iout,*) 'ishell: cmd variable too short.'
write(iout,*) 'Program will stop.'
call flush(iout)
call exit(1)
endif
100 continue
istat=system(cmd)
if(istat.lt.0) then
#else
character*512 cmd
cmd=string // char(0)
istat=system(cmd)
if(istat.ne.0) then
#endif
write(iout,*)
write(iout,*) 'Fatal error in ',string,'.'
write(iout,*) 'Program will stop.'
call flush(iout)
call dmrccend(1)
endif
C
return
end
C
************************************************************************
subroutine dmrccend(ierr)
************************************************************************
* Terminates mrcc
************************************************************************
#include "MRCCCOMMON"
integer ierr
character*16 devnul
common/output/ devnul
#if defined(MPI)
integer mpi_err, mpi_rank
logical mpi_in
#endif
C
write(iout,*)
call system(
$'echo " ************************ "`date +"%F %T"`"' //
$' *************************"'//trim(devnul))
c#if defined(MPI)
c call MPI_Initialized(mpi_in, mpi_err)
c
c if(mpi_in) then
c call MPI_Comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_err)
c
c if(mpi_rank .ne. 0) read(iout, *)
c end if
c#endif
if(ierr.eq.0) then
write(iout,"(22x,a27)") 'Normal termination of mrcc.'
else
write(iout,"(19x,a33)") 'Error at the termination of mrcc.'
endif
write(iout,"(1x,70('*'))")
#if defined(MPI)
call mpifinalize(ierr)
#endif
call delete_file('pids', inbcast)
call exit(ierr)
C
return
end
C
************************************************************************
subroutine getvar(vname,var)
************************************************************************
* Read variables from file VARS
************************************************************************
#include "MRCCCOMMON"
integer i,n
integer*1 var(*)
character*10 vname,cscr
C
open(varsfile,file='VARS',form='UNFORMATTED')
rewind(varsfile)
do
read(varsfile,end=1002) cscr
if(cscr.eq.vname) exit
enddo
backspace(varsfile)
read(varsfile) cscr,n,(var(i),i=1,n)
close(varsfile)
C
return
1002 continue
write(iout,*) 'Corrupted file VARS!'
write(iout,*) 'Variable '//trim(adjustl(vname))//' is missing!'
write(iout,*) 'Please run integ first!'
call mrccend(1)
C
return
end
C
************************************************************************
subroutine getvar_c(vname,var)
************************************************************************
* Read variables from file VARS
************************************************************************
#include "MRCCCOMMON"
integer i,n
character*1 var(*)
character*10 vname,cscr
C
open(varsfile,file='VARS',form='UNFORMATTED')
rewind(varsfile)
do
read(varsfile,end=1002) cscr
if(cscr.eq.vname) exit
enddo
backspace(varsfile)
read(varsfile) cscr,n,(var(i),i=1,n)
close(varsfile)
C
return
1002 continue
write(iout,*) 'Corrupted file VARS!'
write(iout,*) 'Please run integ first!'
call mrccend(1)
C
return
end
C
************************************************************************
subroutine rspmx(h,n,ifile)
************************************************************************
* Read and unpack symmetry-packed matrix
************************************************************************
implicit none
integer n,ifile,i,j,ii
real*8 h(n*n)
C
read(ifile) (h(i),i=1,(n+1)*n/2)
ii=(n+1)*n/2
do j=n,1,-1
do i=n,j,-1
h((j-1)*n+i)=h(ii)
ii=ii-1
enddo
enddo
c call fillup(h,n)
C
return
end
C
************************************************************************
subroutine rspmx2(h,n,ifile)
************************************************************************
* Read symmetry-packed matrix
************************************************************************
implicit none
integer n,ifile,i,j,ii
real*8 h(n*n)
C
read(ifile) (h(i),i=1,(n+1)*n/2)
C
return
end
C
C***********************************************************************
subroutine roeint_kineti(r8heap,i4heap,h,ifile,nbasis)
C***********************************************************************
C Read one-electron quantities
C***********************************************************************
implicit none
integer ii,jj,nn,i,j,ifile,nbasis
real*8 r8heap(*),h(nbasis,nbasis)
integer*4 i4heap(*)
C
call dfillzero(h,nbasis**2)
read(ifile) nn
rewind(ifile)
c write(6,*) nn
read(ifile) ii,(r8heap(jj),jj=1,2*nn)
do ii=1,nn
i=i4heap((ii-1)*4+3)
j=i4heap((ii-1)*4+4)
h(i,j)=r8heap(ii*2-1)
h(j,i)=r8heap(ii*2-1)
c write(6,"(i7,2x,2i4,f12.5)") ii,i,j,h(i,j)
enddo
C
return
end
C
C***********************************************************************
subroutine roeint(r8heap,i4heap,h,ifile,nbasis)
C***********************************************************************
C Read one-electron quantities
C***********************************************************************
implicit none
integer ii,jj,nn,i,j,ifile,nbasis
real*8 r8heap(*),h(nbasis,nbasis)
integer*4 i4heap(*)
C
call dfillzero(h,nbasis**2)
read(ifile) nn
backspace(ifile)
c write(6,*) nn
read(ifile) ii,(r8heap(jj),jj=1,2*nn)
do ii=1,nn
i=i4heap((ii-1)*4+3)
j=i4heap((ii-1)*4+4)
h(i,j)=r8heap(ii*2-1)
h(j,i)=r8heap(ii*2-1)
c write(6,"(i7,2x,2i4,f12.5)") ii,i,j,h(i,j)
enddo
C
return
end
C
C***********************************************************************
subroutine woeint(r8heap,i4heap,h,ifile,itol,nbasis)
C***********************************************************************
C Write one-electron quantities
C***********************************************************************
implicit none
integer ii,jj,nn,i,j,ifile,nbasis
real*8 r8heap(*),h(nbasis,nbasis),itol
integer*4 i4heap(*)
C
nn=0
do i=1,nbasis
do j=1,i
if(dabs(h(i,j)).gt.itol) then
nn=nn+1
r8heap(nn*2-1)=h(i,j)
i4heap((nn-1)*4+3)=i
i4heap((nn-1)*4+4)=j
endif
enddo
enddo
write(ifile) nn,(r8heap(jj),jj=1,2*nn)
C
return
end
C
C***********************************************************************
subroutine woeintu(r8heap,i4heap,h,ifile,itol,nbasis)
C***********************************************************************
C Write one-electron quantities
C***********************************************************************
implicit none
integer ii,jj,nn,i,j,ifile,nbasis
real*8 r8heap(*),h(nbasis,nbasis),itol
integer*4 i4heap(*)
C
nn=0
do i=1,nbasis
do j=1,i
if(dabs(h(j,i)).gt.itol) then
nn=nn+1
r8heap(nn*2-1)=h(j,i)
i4heap((nn-1)*4+3)=j
i4heap((nn-1)*4+4)=i
endif
enddo
enddo
write(ifile) nn,(r8heap(jj),jj=1,2*nn)
C
return
end
C
C***********************************************************************
subroutine wtdmx(r8heap,i4heap,h,ifile,itol,nb1,nb2)
C***********************************************************************
C Write two-dimensional matrix
C***********************************************************************
implicit none
integer ii,jj,nn,i,j,ifile,nb1,nb2
real*8 r8heap(*),h(nb1,nb2),itol
integer*4 i4heap(*)
C
c write(6,"(4f10.5)") h
nn=0
do i=1,nb1
do j=1,nb2
if(dabs(h(i,j)).gt.itol) then
nn=nn+1
r8heap(nn*2-1)=h(i,j)
i4heap((nn-1)*4+3)=i
i4heap((nn-1)*4+4)=j
endif
enddo
enddo
write(ifile) nn,(r8heap(jj),jj=1,2*nn)
C
return
end
C
C***********************************************************************
subroutine wrtmo(r8heap,i4heap,c,ifile,itol,nbf,nbasis)
C***********************************************************************
C Write MO coefficients
C***********************************************************************
implicit none
integer ii,jj,nn,i,j,ifile,nbf,nbasis
real*8 r8heap(*),c(nbf,nbasis),itol
integer*4 i4heap(*)
C
nn=0
do j=1,nbasis
do i=1,nbf
if(dabs(c(i,j)).gt.itol) then
nn=nn+1
r8heap(nn*2-1)=c(i,j)
i4heap((nn-1)*4+3)=i
i4heap((nn-1)*4+4)=j
endif
enddo
enddo
write(ifile) nn,(r8heap(jj),jj=1,2*nn)
C
return
end
C
C***********************************************************************
subroutine rtdmx(r8heap,i4heap,h,ifile,nb1,nb2)
C***********************************************************************
C Read one-electron quantities
C***********************************************************************
implicit none
integer ii,jj,nn,i,j,ifile,nb1,nb2
real*8 r8heap(*),h(nb1,nb2)
integer*4 i4heap(*)
C
call dfillzero(h,nb1*nb2)
read(ifile) nn,(r8heap(jj),jj=1,2*nn)
do ii=1,nn
i=i4heap((ii-1)*4+3)
j=i4heap((ii-1)*4+4)
if(j.le.nb2) h(i,j)=r8heap(ii*2-1)
C write(6,"(2i4,f12.5)") i,j,h(i,j)
enddo
C
return
end
C***********************************************************************
subroutine roequa(r8heap,i4heap,h,ifile,nbasis)
C***********************************************************************
C Read non-symmetric one-electron quantities
C***********************************************************************
implicit none
integer ii,jj,nn,i,j,ifile,nbasis
real*8 r8heap(*),h(nbasis,nbasis)
integer*4 i4heap(*)
C
call dfillzero(h,nbasis**2)
read(ifile) nn,(r8heap(jj),jj=1,2*nn)
do ii=1,nn
i=i4heap((ii-1)*4+3)
j=i4heap((ii-1)*4+4)
h(i,j)=r8heap(ii*2-1)
c write(6,"(2i4,f12.5)") i,j,h(i,j)
enddo
C
return
end
C
C***********************************************************************
subroutine syminv(mat,invmat,n,nn,scr)
C***********************************************************************
C Calculate the inverse of a symmetric real matrix
C The lower triangle of the matrix is used!
C***********************************************************************
implicit none
integer n,nn,i,j
real*8 mat(nn,nn),invmat(nn,nn),scr(n*(n+1)/2+2*n)
C
if (n.gt.1) then
do j=1,n
do i=1,j
scr(i+(j-1)*j/2)=mat(j,i)
enddo
enddo
call dsptrf('U',n,scr,scr(n*(n+1)/2+1),scr(n*(n+1)/2+n+1),i)
call dsptri('U',n,scr,scr(n*(n+1)/2+1),scr(n*(n+1)/2+n+1),i)
if(i.ne.0) then
write(6,*) 'Fatal error at matrix inversion!'
if(i.lt.0) then
write(6,*) 'Illegal arguments!'
else
write(6,*) 'The matrix is singular!'
endif
call mrccend(1)
endif
c write(6,"(10000f9.5)") scr
do j=1,n
do i=1,j
invmat(i,j)=scr(i+(j-1)*j/2)
invmat(j,i)=invmat(i,j)
enddo
enddo
elseif (n.eq.1) then !NP
if(dabs(mat(1,1)).lt.1.d-15) then
write(6,*) 'The matrix for inversion has 1 element:',mat(1,1)
call mrccend(1)
endif
invmat(1,1)=1.d0/mat(1,1)
endif
C
return
end
C
C***********************************************************************
subroutine syminvpd(mat,n,nn)
C***********************************************************************
C Calculate the inverse of a symmetric real positive definite matrix
C The lower triangle of the matrix is used!
C***********************************************************************
implicit none
integer n,nn,i
integer*4 ii
equivalence(i,ii)
real*8 mat(nn,nn)
C
call dpotrf('L',n,mat,nn,ii)
if(ii.ne.0) goto 9876
call dpotri('L',n,mat,nn,ii)
if(ii.ne.0) goto 9876
C
return
9876 write(6,*) 'Fatal error at matrix inversion!'
write(6,*) 'The matrix is not positive definite!'
call mrccend(1)
end
C
************************************************************************
subroutine invsqrt(mat,n,scr,iout,erri,tol,iflag)
************************************************************************
* Calculate the inverse/square root of a symmetric, positive def. matrix
* The lower triangle of the matrix is stored
* iflag=1 - inverse square root
* iflag=2 - square root
* iflag=3 - Cholesky deco of the inverse
* iflag=4 - inverse
************************************************************************
implicit none
integer n,iout,i,j,erri,iflag,ii,info8
integer*4 err,info4
real*8 mat(n,n),scr(*),ss,tol,mr
equivalence(info8,info4)
integer imem,imem1,maxcor,memfree,memreq
common/actmem/ imem,maxcor
common/memcom/ imem1
C
memfree = maxcor - (imem - imem1)
c write(6,*)
c do i=1,n
c write(6,"(10000f9.5)") (mat(i,j),j=1,n)
c enddo
if(iflag.ge.3) then
C Inverse of the matrix
call dpotrf('L',n,mat,n,info4)
if(info4.ne.0) then
write(iout,*) 'Fatal error at Cholesky decomposition!'
call mrccend(1)
endif
call dpotri('L',n,mat,n,info4)
if(info4.ne.0) then
write(iout,*) 'Fatal error at matrix inversion!'
call mrccend(1)
endif
if(iflag.eq.4) then
erri=0
return
endif
C Cholesky deco of the inverse of the matrix
call dpotrf('L',n,mat,n,info4)
if(info4.ne.0) then
write(iout,*) 'Fatal error at Cholesky decomposition!'
call mrccend(1)
endif
do j=2,n
mat(1:j-1,j)=0.d0
enddo
erri=0
c call dtrtri('L','N',n,mat,n,i)
c if(i.ne.0) then
c write(iout,*) 'Fatal error at matrix inversion!'
c call mrccend(1)
c endif
c do i=1,n
c do j=1,min(i,n-i)
c ss=mat(i,j)
c mat(i,j)=mat(n-j+1,n-i+1)
c mat(n-j+1,n-i+1)=ss
c enddo
c enddo
c write(6,*)
c do i=1,n
c write(6,"(10000f9.5)") (mat(i,j),j=1,n)
c enddo
c return
C
c ii=0
c do i=1,n
c do j=1,i
c ii=ii+1
c scr(ii)=mat(i,j)
c enddo
c enddo
c call dsptrf('U',n,scr,scr(n*(n+1)/2+1),scr(n*(n+1)/2+n+1),i)
c call dsptri('U',n,scr,scr(n*(n+1)/2+1),scr(n*(n+1)/2+n+1),i)
c if(i.ne.0) then
c write(iout,*) 'Fatal error at matrix inversion!'
c call mrccend(1)
c endif
c call dpptrf('U',n,scr,i)
c if(i.ne.0) then
c write(iout,*) 'Fatal error at Cholesky decomposition!'
c call mrccend(1)
c endif
c ii=0
c do i=1,n
c do j=1,i
c ii=ii+1
c mat(i,j)=scr(ii)
c enddo
c do j=i+1,n
c mat(i,j)=0.d0
c enddo
c enddo
c erri=0
c write(6,*)
c do i=1,n
c write(6,"(10000f9.5)") (mat(i,j),j=1,n)
c enddo
C dpotf2.f dpptrf.f dpotrf.f
else
C Inverse/square root
call dsyev('V','L',n,mat,n,scr,mr,-1,err)
c memreq = idnint(mr)
memreq = int(mr)
if(memfree .le. memreq) then
write(iout, '(a)') ' Error: Insufficient memory'
write(iout, '(i8,a)') (memreq-memfree)*8/1024**2,
$ ' MB more memory needed'
call mrccend(1)
end if
c write(6,*) (mat(i,i),i=1,n)
c write(6,"(9f10.6)") mat
C Diagonalize
call dsyev('V','L',n,mat,n,scr,scr(n+1),3*n**2,err)
c write(6,*) 'eig'
c write(6,"(10000f20.12)") (scr(j),j=1,n)
erri=err
if(err.ne.0) then
write(iout,*) 'Inverse square root: Fatal error at the ' //
$'diagonalization of the matrix!'
call mrccend(1)
endif
C Inverse square root
call dcopy(n**2,mat,1,scr(n+1),1)
do i=1,n
ss=scr(i)
if(ss.lt.tol) then
c if(ss.lt.1d-10) then
write(iout,*)
$'Inverse square root: The matrix is not positive definite!'
c write(iout,"(' Critical eigenvalue: ',1pe8.2)") ss
erri=1
c call mrccend(1)
return
else
c if(iflag.eq.4) then
c ss=1.d0/ss !inverse
c else
ss=dsqrt(ss) !square root
if(iflag.eq.1) ss=1.d0/ss !inverse square root
c endif
call dscal(n,ss,mat(1,i),1)
endif
enddo
call dgemm('n','t',n,n,n,1.d0,mat,n,scr(n+1),n,0.d0,
$scr(n+n**2+1),n)
call dcopy(n**2,scr(n+n**2+1),1,mat,1)
c write(6,"(2f10.5)") mat
endif
C
return
end
C
************************************************************************
subroutine sppack(n,m,c,spc,spi,tol)
************************************************************************
* Pack sparse matrix
************************************************************************
implicit none
integer n,m,i,j,nn
integer*8 spi(0:3*n*m)
real*8 c(n,m),spc(0:3*n*m),tol
C
nn=0
do i=1,n
c write(6,"(10000f14.8)") (c(i,j),j=1,m)
do j=1,m
if(dabs(c(i,j)).gt.tol) then
spc(3*nn+1)=c(i,j)
spi(3*nn+2)=i
spi(3*nn+3)=j
nn=nn+1
endif
enddo
enddo
spi(0)=3*nn-3
C
return
end
C
************************************************************************
subroutine spgemmn(m,n,k,spa,spi,b,c)
************************************************************************
* Multiplication of sparse matrix A with matrix B
************************************************************************
implicit none
integer m,n,i,j,k,l,nn
integer*8 spi(0:3*m*k)
real*8 spa(0:3*m*k),b(k,n),c(m,n),sum
C
c call dfillzero(c,m*n)
c=0.d0
do nn=0,spi(0),3
sum=spa(nn+1)
i= spi(nn+2)
l= spi(nn+3)
do j=1,n
c(i,j)=c(i,j)+sum*b(l,j)
enddo
enddo
C
return
end
C
************************************************************************
subroutine spgemmt(m,n,k,spa,spi,b,c)
************************************************************************
* Multiplication of sparse matrix A with matrix B^T
************************************************************************
implicit none
integer m,n,i,j,k,l,nn
integer*8 spi(0:3*m*k)
real*8 spa(0:3*m*k),b(n,k),c(m,n),sum
C
c call dfillzero(c,m*n)
c=0.d0
do nn=0,spi(0),3
sum=spa(nn+1)
i= spi(nn+2)
l= spi(nn+3)
c c(i,1:n)=c(i,1:n)+sum*b(1:n,l)
do j=1,n
c(i,j)=c(i,j)+sum*b(j,l)
enddo
enddo
C
return
end
C
************************************************************************
subroutine spgemmr(m,n,k,spa,spi,b,c)
************************************************************************
* Multiplication of matrix B with the transpose of sparse matrix A
************************************************************************
implicit none
integer m,n,i,j,k,l,nn
integer*8 spi(0:3*n*k)
real*8 spa(0:3*n*k),b(m,k),c(m,n),sum
C
c call dfillzero(c,m*n)
c=0.d0
do nn=0,spi(0),3
sum=spa(nn+1)
j= spi(nn+2)
l= spi(nn+3)
do i=1,m
c(i,j)=c(i,j)+b(i,l)*sum
enddo
enddo
C
return
end
C
************************************************************************
subroutine mxtrans(b,n,m,file,irec,inewadd)
************************************************************************
* General matrix transposition, the matrix is read from file *
************************************************************************
#include "MRCCCOMMON"
integer n,m,i,j,ii,file,irec,inewadd,ifrec
real*8 b(m,n)
C
if(m*n.le.0) return
ifrec=irec+(inewadd-mod(inewadd,ibufln))/ibufln
read(file,rec=ifrec) ibuf
ii=mod(inewadd,ibufln)
do j=1,m
do i=1,n
ii=ii+1
b(j,i)=ibuf(ii)
if(ii.eq.ibufln) then
ii=0
ifrec=ifrec+1
if(i*j.lt.n*m) read(file,rec=ifrec) ibuf
endif
enddo
enddo
C
return
end
C
************************************************************************
subroutine gtrans(a,b,n,m)
************************************************************************
* General matrix transposition
************************************************************************
implicit none
integer n,m,i,j
real*8 a(n,m),b(m,n)
C
do i=1,n
do j=1,m
b(j,i)=a(i,j)
enddo
enddo
C
return
end
C
************************************************************************
subroutine fillup(a,n)
************************************************************************
* Fill up the upper triangle of a symmetric matrix
************************************************************************
implicit none
integer n,i,j
real*8 a(n,n)
C
do i=1,n
do j=1,i-1
a(j,i)=a(i,j)
enddo
enddo
C
return
end
C
************************************************************************
subroutine filllo(a,n)
************************************************************************
* Fill up the lower triangle of a symmetric matrix
************************************************************************
implicit none
integer n,i,j
real*8 a(n,n)
C
do i=1,n
do j=1,i-1
a(i,j)=a(j,i)
enddo
enddo
C
return
end
C
************************************************************************
subroutine filezero(file,ilo,iup)
************************************************************************
* Initialize files *
************************************************************************
#include "MRCCCOMMON"
integer file,ilo,iup,i
C
call dfillzero(ibuf,ibufln)
do i=ilo,iup
write(file,rec=i) ibuf
enddo
C
return
end
C
C***********************************************************************
subroutine readmo(r8heap,i4heap,mo,ifile,nbf,nbasis)
C***********************************************************************
C Read MO coefficients (without core orbitals)
C***********************************************************************
implicit none
integer ii,jj,nn,i,j,ifile,nbasis,nbf,ncore
real*8 r8heap(*),mo(nbf,nbasis)
integer*4 i4heap(*)
C
call dfillzero(mo,nbf*nbasis)
ncore=nbf-nbasis
read(ifile) nn,(r8heap(jj),jj=1,2*nn)
do ii=1,nn
i=i4heap((ii-1)*4+3)
j=i4heap((ii-1)*4+4)
if(j.gt.ncore) mo(i,j-ncore)=r8heap(ii*2-1)
enddo
C
return
end
C
************************************************************************
subroutine gbcheck(lgb,atsymbas,gbasfile)
************************************************************************
* Checks if file GENBAS exists and contains a basis set/grid
************************************************************************
implicit none
integer gbasfile
character*24 atsymbas,glabel
logical lgb
C
inquire(file='GENBAS',exist=lgb)
if(lgb) then
open(gbasfile,file='GENBAS')
glabel=' '
do while(trim(atsymbas).ne.trim(glabel))
read(gbasfile,*,end=1000) glabel
call lowercase(glabel,glabel,24)
enddo
c close(gbasfile)
endif
C
return
1000 continue
close(gbasfile)
lgb=.false.
end
C
************************************************************************
subroutine motransp(nmatrix,nbasis,mo1,mo2,lll)
************************************************************************
* Transpose MO coefficient matrix
************************************************************************
implicit none
integer nmatrix,nbasis,i,mu
real*8 mo1(nmatrix,nbasis),mo2(nbasis,nmatrix)
logical lll
C
if(lll) then
do i=1,nmatrix
do mu=1,nbasis
mo2(mu,i)=mo1(i,mu)
enddo
enddo
else
do i=1,nmatrix
do mu=1,nbasis
mo1(i,mu)=mo2(mu,i)
enddo
enddo
endif
C
return
end
C
************************************************************************
SUBROUTINE QSORTD (ORD,N,A)
************************************************************************
* SORTS THE ARRAY A(I),I=1,2,...,N BY PUTTING THE
* ASCENDING ORDER VECTOR IN ORD. THAT IS ASCENDING ORDERED A
* IS A(ORD(I)),I=1,2,...,N; DESCENDING ORDER A IS A(ORD(N-I+1)),
* I=1,2,...,N . THIS SORT RUNS IN TIME PROPORTIONAL TO N LOG N .
************************************************************************
IMPLICIT INTEGER (A-Z)
C
DIMENSION ORD(N),POPLST(2,20)
REAL*8 X,XX,Z,ZZ,Y
C
C TO SORT DIFFERENT INPUT TYPES, CHANGE THE FOLLOWING
C SPECIFICATION STATEMENTS; FOR EXAMPLE, FOR FORTRAN CHARACTER
C USE THE FOLLOWING: CHARACTER *(*) A(N)
C
REAL*8 A(N)
C
NDEEP=0
U1=N
L1=1
DO 1 I=1,N
1 ORD(I)=I
2 IF (U1.LE.L1) RETURN
C
3 L=L1
U=U1
C
C PART
C
4 P=L
Q=U
C FOR CHARACTER SORTS, THE FOLLOWING 3 STATEMENTS WOULD BECOME
C X = ORD(P)
C Z = ORD(Q)
C IF (A(X) .LE. A(Z)) GO TO 2
C
C WHERE "CLE" IS A LOGICAL FUNCTION WHICH RETURNS "TRUE" IF THE
C FIRST ARGUMENT IS LESS THAN OR EQUAL TO THE SECOND, BASED ON "LEN"
C CHARACTERS.
C
X=A(ORD(P))
Z=A(ORD(Q))
IF (X.LE.Z) GO TO 5
Y=X
X=Z
Z=Y
YP=ORD(P)
ORD(P)=ORD(Q)
ORD(Q)=YP
5 IF (U-L.LE.1) GO TO 15
XX=X
IX=P
ZZ=Z
IZ=Q
C
C LEFT
C
6 P=P+1
IF (P.GE.Q) GO TO 7
X=A(ORD(P))
IF (X.GE.XX) GO TO 8
GO TO 6
7 P=Q-1
GO TO 13
C
C RIGHT
C
8 Q=Q-1
IF (Q.LE.P) GO TO 9
Z=A(ORD(Q))
IF (Z.LE.ZZ) GO TO 10
GO TO 8
9 Q=P
P=P-1
Z=X
X=A(ORD(P))
C
C DIST
C
10 IF (X.LE.Z) GO TO 11
Y=X
X=Z
Z=Y
IP=ORD(P)
ORD(P)=ORD(Q)
ORD(Q)=IP
11 IF (X.LE.XX) GO TO 12
XX=X
IX=P
12 IF (Z.GE.ZZ) GO TO 6
ZZ=Z
IZ=Q
GO TO 6
C
C OUT
C
13 CONTINUE
IF (.NOT.(P.NE.IX.AND.X.NE.XX)) GO TO 14
IP=ORD(P)
ORD(P)=ORD(IX)
ORD(IX)=IP
14 CONTINUE
IF (.NOT.(Q.NE.IZ.AND.Z.NE.ZZ)) GO TO 15
IQ=ORD(Q)
ORD(Q)=ORD(IZ)
ORD(IZ)=IQ
15 CONTINUE
IF (U-Q.LE.P-L) GO TO 16
L1=L
U1=P-1
L=Q+1
GO TO 17
16 U1=U
L1=Q+1
U=P-1
17 CONTINUE
IF (U1.LE.L1) GO TO 18
C
C START RECURSIVE CALL
C
NDEEP=NDEEP+1
POPLST(1,NDEEP)=U
POPLST(2,NDEEP)=L
GO TO 3
18 IF (U.GT.L) GO TO 4
C
C POP BACK UP IN THE RECURSION LIST
C
IF (NDEEP.EQ.0) GO TO 2
U=POPLST(1,NDEEP)
L=POPLST(2,NDEEP)
NDEEP=NDEEP-1
GO TO 18
C
C END SORT
C END QSORT
C
END
C
!NP
c {{{ subroutines timeadd time0
subroutine timeadd(t,t0)
implicit none
real*8 t(6),tmp,t0(2)
integer w,ps
c integer omp_get_thread_num,
c $ omp_get_ancestor_thread_num,omp_get_level
c write(*,*) "timeadd", omp_get_thread_num(),
c $ omp_get_ancestor_thread_num(omp_get_level()-1)
c$OMP CRITICAL
call CPU_TIME(tmp)
call SYSTEM_CLOCK(w,ps)
t(1)=tmp-t0(1)
t(2)=dble(w)/dble(ps)-t0(2)
c sum of times
t(3)=t(3)+t(1)
t(4)=t(4)+t(2)
c max times
t(5)=max(t(1),t(5))
t(6)=max(t(2),t(6))
c$OMP END CRITICAL
return
end
subroutine time0(t)
implicit none
real*8 t(2)
integer w,ps
call CPU_TIME(t(1))
call SYSTEM_CLOCK(w,ps)
t(2)=dble(w)/dble(ps)
return
end
c }}}
!NP
************************************************************************
subroutine dsymv_sp(n,a,x,y)
************************************************************************
* y=Ax matrix-vector product, symmetry-packed matrix, lower triangle
************************************************************************
implicit none
integer n,i,j,ii,nn
real*8 a((n+1)*n/2),x(n),y(n),tmp1,tmp2,tmp
c real*8 a(n,n),x(n),y(n)
C
y=0.d0
nn=2*n+2
do j=1,n
tmp1=x(j)
tmp2=0.d0
ii=(nn-j)*(j-1)/2-j+1
y(j)=y(j)+tmp1*a(ii+j)
C$OMP PARALLEL DO Schedule(Static)
C$OMP& DEFAULT(SHARED)
C$OMP& PRIVATE(i,tmp)
C$OMP& REDUCTION(+:tmp2)
do i=j+1,n
tmp=a(ii+i)
y(i)=y(i)+tmp1*tmp
tmp2=tmp2+tmp*x(i)
enddo
C$OMP END PARALLEL DO
y(j)=y(j)+tmp2
enddo
C
return
end
C
************************************************************************
subroutine getenergy(ener,cscr15)
************************************************************************
* Get energy from iface file
************************************************************************
#include "MRCCCOMMON"
integer i,j
real*8 ener
character(len=8) cscr8
character(len=15) cscr15
C
open(unit=ifcfile,file='iface',status='old')
i=0
do
read(ifcfile,*,end=7597)
i=i+1
enddo
7597 rewind(ifcfile)
do j=1,i-1
read(ifcfile,*)
enddo
read(ifcfile,7596) cscr8,cscr15,i,i,i,ener
close(ifcfile)
7596 format(a8,1x,a15,1x,3(i2,1x),1pe23.15)
C
return
end
C
************************************************************************
logical function lanycc(calc)
************************************************************************
* Check if calc refers to any of the CC or CI methods computed with mrcc
* as well as ADC(2) and CIS(D_Inf)
* Warning: MP2, spin scaled MP2 and all RPA variants are left out intentionally!
************************************************************************
implicit none
character*16 calc,xcalc
character*1 xcalc1(16)
equivalence(xcalc,xcalc1)
integer i,ii
C
xcalc=calc
do i=2,15
ii=ichar(xcalc1(i))
if(ii.ge.49.and.ii.le.57.and.(xcalc1(i-1).eq.'['.or.
$xcalc1(i-1).eq.'('.or.xcalc1(i+1).eq.']'.or.xcalc1(i+1).eq.')'.or.
$xcalc1(i-1).eq.'c'.or.xcalc1(i-1).eq.'n')) xcalc1(i)='n'
enddo
do i=2,13
if(xcalc1(i ).eq.'-'.and.
$ xcalc1(i+1).eq.'f'.and.
$ xcalc1(i+2).eq.'1'.and.
$ xcalc1(i+3).eq.'2') then
xcalc1(i:i+3)=' '
if(xcalc1(1).eq.'c'.and.xcalc1(2).eq.'i') then
write(6,*) 'Explicit correlation is not available for CI!'
call mrccend(1)
endif
exit
endif
enddo
c
lanycc=
$ xcalc.eq.'mp3 '.or.
$ xcalc.eq.'ccs '.or.
$ xcalc.eq.'ccsd '.or.
$ xcalc.eq.'ccsdt '.or.
$ xcalc.eq.'ccsdtq '.or.
$ xcalc.eq.'ccsdtqp '.or.
$ xcalc.eq.'cc(n) '.or.
$ xcalc.eq.'cc(nn) '.or.
$ xcalc.eq.'ccsd[t] '.or.
$ xcalc.eq.'ccsdt[q] '.or.
$ xcalc.eq.'ccsdtq[p] '.or.
$ xcalc.eq.'cc(n)[n] '.or.
$ xcalc.eq.'cc(n)[nn] '.or.
$ xcalc.eq.'cc(nn)[nn] '.or.
$ xcalc.eq.'ccsd(t) '.or.
$ xcalc.eq.'ccsdt(q) '.or.
$ xcalc.eq.'ccsdtq(p) '.or.
$ xcalc.eq.'cc(n)(n) '.or.
$ xcalc.eq.'cc(n)(nn) '.or.
$ xcalc.eq.'cc(nn)(nn) '.or.
$ xcalc.eq.'ccsd(t)_l '.or.
$ xcalc.eq.'ccsdt(q)_l '.or.
$ xcalc.eq.'ccsdtq(p)_l '.or.
$ xcalc.eq.'cc(n)(n)_l '.or.
$ xcalc.eq.'cc(n)(nn)_l '.or.
$ xcalc.eq.'cc(nn)(nn)_l '.or.
$ xcalc.eq.'ccsdt-1a '.or.
$ xcalc.eq.'ccsdtq-1a '.or.
$ xcalc.eq.'ccsdtqp-1a '.or.
$ xcalc.eq.'cc(n)-1a '.or.
$ xcalc.eq.'cc(nn)-1a '.or.
$ xcalc.eq.'ccsdt-1b '.or.
$ xcalc.eq.'ccsdtq-1b '.or.
$ xcalc.eq.'ccsdtqp-1b '.or.
$ xcalc.eq.'cc(n)-1b '.or.
$ xcalc.eq.'cc(nn)-1b '.or.
$ xcalc.eq.'ccn '.or.
$ xcalc.eq.'ccnn '.or.
$ xcalc.eq.'ccsdt-3 '.or.
$ xcalc.eq.'ccsdtq-3 '.or.
$ xcalc.eq.'ccsdtqp-3 '.or.
$ xcalc.eq.'cc(n)-3 '.or.
$ xcalc.eq.'cc(nn)-3 '.or.
$ xcalc.eq.'ccsdt[q]/a '.or.
$ xcalc.eq.'ccsdtq[p]/a '.or.
$ xcalc.eq.'cc(n)[n]/a '.or.
$ xcalc.eq.'cc(n)[nn]/a '.or.
$ xcalc.eq.'cc(nn)[nn]/a '.or.
$ xcalc.eq.'ccsdt[q]/b '.or.
$ xcalc.eq.'ccsdtq[p]/b '.or.
$ xcalc.eq.'cc(n)[n]/b '.or.
$ xcalc.eq.'cc(n)[nn]/b '.or.
$ xcalc.eq.'cc(nn)[nn]/b '.or.
$ xcalc.eq.'ccsdt(q)/a '.or.
$ xcalc.eq.'ccsdtq(p)/a '.or.
$ xcalc.eq.'cc(n)(n)/a '.or.
$ xcalc.eq.'cc(n)(nn)/a '.or.
$ xcalc.eq.'cc(nn)(nn)/a '.or.
$ xcalc.eq.'ccsdt(q)/b '.or.
$ xcalc.eq.'ccsdtq(p)/b '.or.
$ xcalc.eq.'cc(n)(n)/b '.or.
$ xcalc.eq.'cc(n)(nn)/b '.or.
$ xcalc.eq.'cc(nn)(nn)/b '.or.
$ xcalc.eq.'cis '.or.
$ xcalc.eq.'cisd '.or.
$ xcalc.eq.'cisdt '.or.
$ xcalc.eq.'cisdtq '.or.
$ xcalc.eq.'cisdtqp '.or.
$ xcalc.eq.'fci '.or.
$ xcalc.eq.'ci(n) '.or.
$ xcalc.eq.'adc(n) '.or. ! or calc.eq.'adc(2) ' but wrong with xcalc.eq.'adc(2) '
$ xcalc.eq.'cis(di) '.or.
$ xcalc.eq.'ci(nn) '
return
end
C
c {{{ subroutine setenvcalc
subroutine setenvcalc(envcalc,envcalcfull,minpfile,iout,lrc)
************************************************************************
* set calc keyword for the environment
************************************************************************
implicit none
integer clen
parameter(clen=16)
character*1 c16_1(clen)
character*16 envcalc,c16,envcalcfull
character*32 c32
character*8 cscr8
integer minpfile,i,iout
equivalence(c16,c16_1)
logical localmp2,lrc,lstartwlno,lstartwl
c
open(minpfile,file='MINP')
call getkeym('corembed',8,cscr8,8)
read(minpfile,*)
read(minpfile,*) c16
close(minpfile)
call lowercase(c16,c16,clen)
c32(1:16)=c16
c32(17:32)=" "
lrc=.true.
if (lstartwlno(c32)) then ! for the LNO- synonym
call removeprefix(c16_1,clen,4)
elseif (lstartwl(c32)) then ! for the L synonym
call removeprefix(c16_1,clen,1)
elseif (trim(c16).eq.'hf') then
c16='scf '
lrc=.false.
elseif (trim(c16).eq.'hf+lrc') then
c16='scf '
else
write(iout,*)'Input error: only local correlation methods and
$HF or HF+lrc are available as environment models of corembed'
call mrccend(1)
endif
envcalcfull=c16
localmp2=(envcalcfull.eq.'mp2 '.or.
$ envcalcfull.eq.'scs-mp2 '.or.
$ envcalcfull.eq.'sos-mp2 ')
envcalc=envcalcfull
if (localmp2) envcalc='mp2 '
c
return
end
c }}}
C
c {{{ function setccprog
character*4 function setccprog(envcalc)
************************************************************************
* set ccprog keyword for the environment
************************************************************************
implicit none
character*16 envcalc
logical domrcc,lanycc
c
domrcc=.true.
if (envcalc.eq.'ccsd '.or.
$ envcalc.eq.'ccsd(t) ') then
setccprog='ccsd'
domrcc=.false.
endif
if (lanycc(envcalc).and.domrcc) then
setccprog='mrcc'
endif
c
return
end
c }}}
C
************************************************************************
real*8 function expval(o,d,n)
************************************************************************
* Calculate the expectation value of a one-electron operator
************************************************************************
implicit none
integer n,i,j
real*8 o(n,n),d(n,n),ss
C
ss=0.d0
C$OMP PARALLEL DO Schedule(Dynamic)
C$OMP& DEFAULT(SHARED)
C$OMP& PRIVATE(j)
C$OMP& REDUCTION(+:ss)
do j=1,n
ss=ss+o(j,j)*d(j,j)+2.d0*dot_product(o(j+1:n,j),d(j+1:n,j))
enddo
C$OMP END PARALLEL DO
expval=ss
C
return
end
C
************************************************************************
subroutine getgrfile(grfile,atsymbol,basname,charge,ich,natoms,
$scfiguess)
************************************************************************
* Construct the name of grfile
************************************************************************
implicit none
integer charge,ich,natoms,i
character*1 bname1(20)
character*2 ch2,atsymbol
character*7 scfiguess
character*16 cscr16
character*20 basname,bname
character*40 grfile,c40
equivalence(bname,bname1)
C
if(mod(charge,natoms).eq.0) then
ich=charge/natoms
else
ich=0
endif
write(c40,'(i5)') ich
bname=basname
do i=1,20
if(bname1(i).eq.'('.or.bname1(i).eq.')') bname1(i)='-'
if(bname1(i).eq.'*') bname1(i)='x'
enddo
if(scfiguess.eq.'gfnrest') then
grfile='SCFCOREDENSITIES.' // trim(adjustl(atsymbol)) //
$'.'// trim(adjustl(bname)) // '.' // trim(adjustl(c40))
else !if(scfiguess.eq.'sad') then
grfile='SCFDENSITIES.' // trim(adjustl(atsymbol)) //
$'.'// trim(adjustl(bname)) // '.' // trim(adjustl(c40))
endif
C
return
end
C
#if defined (LIBXC)
************************************************************************
subroutine libxcconv(dft)
************************************************************************
* Convert mrcc functional name to Libxc identifier
************************************************************************
implicit none
character(len=32) dft
C
select case(trim(dft))
C LDA
case('pz')
dft='LDA_C_PZ'
case('vwn1')
dft='LDA_C_VWN_1'
case('vwn2')
dft='LDA_C_VWN_2'
case('vwn3')
dft='LDA_C_VWN_RPA'
case('vwn4')
dft='LDA_C_VWN_4'
C GGA
case('g96')
dft='GGA_X_G96'
case('mpw91x')
dft='GGA_X_MPW91'
case('pbeh')
dft='GGA_X_WPBEH'
c case('vv10')
c dft='GGA_XC_VV10'
case('mpwlyp1w')
dft='GGA_XC_MPWLYP1W'
case('xlyp')
dft='GGA_XC_XLYP'
C Hybrid GGA
case('x3lyp')
dft='HYB_GGA_XC_X3LYP'
case('b1lyp')
dft='HYB_GGA_XC_B1LYP'
case('o3lyp')
dft='HYB_GGA_XC_O3LYP'
C Meta-GGA
case('b95')
dft='MGGA_C_BC95'
case('revtpssx')
dft='MGGA_X_REVTPSS'
case('revtpssc')
dft='MGGA_C_REVTPSS'
case('tpssx')
dft='MGGA_X_TPSS'
case('tpssc')
dft='MGGA_C_TPSS'
case('scanx')
dft='MGGA_X_SCAN'
case('scanc')
dft='MGGA_C_SCAN'
case('revscanc')
dft='MGGA_C_REVSCAN'
C Hybrid meta-GGA
case('tpssh')
dft='HYB_MGGA_XC_TPSSH'
case('revtpssh')
dft='HYB_MGGA_XC_REVTPSSH'
case('mpw1b95')
dft='HYB_MGGA_XC_MPW1B95'
case('pw6b95')
dft='HYB_MGGA_XC_PW6B95'
end select
call lowercase(dft,dft,32)
C
return
end
C
#endif
************************************************************************
subroutine dmatmul(transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,
$ldc,choice,nthread)
************************************************************************
* real*8 matrix multiplication with intrinsic matmul or dgemm
* see Lapack dgemm for the explanation of variables
* beta.ne.0.d0 is not implemented for matmul
************************************************************************
implicit none
character*6 choice,finalc
character*1 transa,transb
integer m,n,k,lda,ldb,ldc,stacksize,nthread
parameter (stacksize=6*10**5)
real*8 alpha,beta,a(lda,*),b(ldb,*),c(ldc,*)
c real*8 alpha,beta,a(m,k),b(k,n),c(m,n)
c
finalc='matmul'
if (choice.eq.'auto '.or.choice.eq.'matmul') then
c stack owerflow occurs with ifort + matmul for large matrices
if ((m*k+n*k+m*n)*nthread.gt.stacksize.or.beta.ne.0.d0)
$ finalc='dgemm '
endif
c
if (choice.eq.'dgemm '.or.finalc.eq.'dgemm ') then
call dgemm(transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc)
else
if (transa.eq.'n'.and.transb.eq.'n') then
c(1:m,1:n)=alpha*matmul(a(1:m,1:k),b(1:k,1:n))
elseif (transa.eq.'t'.and.transb.eq.'n') then
c(1:m,1:n)=alpha*matmul(transpose(a(1:k,1:m)),b(1:k,1:n))
elseif (transa.eq.'n'.and.transb.eq.'t') then
c(1:m,1:n)=alpha*matmul(a(1:m,1:k),transpose(b(1:n,1:k)))
elseif (transa.eq.'t'.and.transb.eq.'t') then
c(1:m,1:n)=alpha*
$ matmul(transpose(a(1:k,1:m)),transpose(b(1:n,1:k)))
endif
endif
c
return
end
C
************************************************************************
subroutine symmat(a,n)
************************************************************************
* Symmetrize nxn matrix a
************************************************************************
implicit none
integer n,i,j
real*8 a(n,n),tmp
C
do i=1,n
do j=1,i-1
tmp=0.5d0*(a(i,j)+a(j,i))
a(i,j)=tmp
a(j,i)=tmp
enddo
enddo
C
return
end
C
************************************************************************
subroutine symmat2(a,n)
************************************************************************
* Symmetrize nxn matrix a without the 0.5 factor
************************************************************************
implicit none
integer n,i,j
real*8 a(n,n),tmp
C
do i=1,n
do j=1,i
tmp=a(i,j)+a(j,i)
a(i,j)=tmp
a(j,i)=tmp
enddo
enddo
C
return
end
C
************************************************************************
subroutine gtrmat(trmat,nispher,nicontr,nicart,niprim,nprimmax,
$ncontrmax,gcoef,ctostr)
************************************************************************
* Genarate Cartesian-primitive -> spherical-contracted transformation
************************************************************************
implicit none
integer nispher,nicontr,nicart,niprim,nprimmax,ncontrmax
integer ispher,icontr,icart,iprim,ii
real*8 gcoef(nprimmax,ncontrmax)
real*8 ctostr(nispher,nicart)
real*8 trmat(nispher*nicontr*nicart*niprim)
C
ii=0
do iprim=1,niprim
do icart=1,nicart
do icontr=1,nicontr
do ispher=1,nispher
ii=ii+1
trmat(ii)=ctostr(ispher,icart)*gcoef(iprim,icontr)
enddo
enddo
enddo
enddo
C
return
end
C
************************************************************************
character*1 function angsym(l)
************************************************************************
* Genarate symbol for azimuthal quantum number
************************************************************************
implicit none
integer l
C
if(l.eq.0) then
angsym='s'
else if(l.eq.1) then
angsym='p'
else if(l.eq.2) then
angsym='d'
else if(l.eq.3) then
angsym='f'
else if(l.eq.4) then
angsym='g'
else if(l.eq.5) then
angsym='h'
else if(l.eq.6) then
angsym='i'
else
angsym=char(l+100)
endif
C
return
end
c
************************************************************************
subroutine moldenloc(mocoeffull,eo,ev,eob,evb,perm,nbasis,
$ moldenfile,nbf,no,nv,nob,nvb,eigenvalue,scrfile1,oslcc)
************************************************************************
* Print (strictly) local MO's (e.g. in an extended domain) or NO's in cis.f
************************************************************************
C
implicit none
C
integer nbasis,i,j,ii,jj,nn,perm(nbf),mocoeffile,moldenfile,nbf
integer no,nv,nob,nvb,scrfile1
real*8 eigenvalue(nbasis),eo(no),ev(nv),eob(nob),evb(nvb)
real*8 mocoeffull(nbf,nbasis)
character*4 mark
character*15 cscr
logical oslcc
C
do i=1,no
eigenvalue(i)=eo(i)
enddo
do i=1,nv
eigenvalue(no+i)=ev(i)
enddo
if (oslcc) then
do i=1,nob
eigenvalue(no+nv+i) = eob(i)
enddo
do i=1,nvb
eigenvalue(no+nv+nob+i) = evb(i)
enddo
endif
C Write MO coefficients to moldenfile
rewind(moldenfile)
read(moldenfile,"(a15,1000000i6)") cscr,(perm(i),i=1,nbf)
open(scrfile1,file='MOLDEN.perm')
read(scrfile1,*,err=101) perm(1:nbf)
101 continue
close(scrfile1)
mark='aaaa'
do while (mark.ne.'[MO]')
read(moldenfile,"(4a)")mark
enddo
do i=1,no+nv
write(moldenfile,"(' Ene=',f14.4)") min(eigenvalue(i),1.d7)
write(moldenfile,"(12a)") ' Spin= Alpha'
if(i.le.no.and..not.oslcc) then
write(moldenfile,"(' Occup=',f14.4)") 2.d0
elseif ((i.le.no.and.oslcc)) then
write(moldenfile,"(' Occup=',f14.4)") 1.d0
elseif (i.gt.no) then
write(moldenfile,"(' Occup=',f14.4)") 0.d0
endif
do j=1,nbf
write(moldenfile,"(i4,f18.10)") j,mocoeffull(perm(j),i)
enddo
enddo
if (oslcc) then
do i=1,nob+nvb
write(moldenfile,"(' Ene=',f14.4)")
$ min(eigenvalue(no+nv+i),1.d7)
write(moldenfile,"(11a)") ' Spin= Beta'
if(i.le.nob) then
write(moldenfile,"(' Occup=',f14.4)") 1.d0
else
write(moldenfile,"(' Occup=',f14.4)") 0.d0
endif
do j=1,nbf
write(moldenfile,"(i4,f18.10)") j,mocoeffull(perm(j),no+nv+i)
enddo
enddo
endif !oslcc
C
end
C
************************************************************************
subroutine moldeninp(natoms,nangmax,ncontrmax,nprimmax,atsymbol,
$nang,ncontr,nprim,gexp,gcoef,coord,atnum,cartg,molperm)
************************************************************************
* Write MOLDEN input file: coordinates and AO basis info
************************************************************************
#include "MRCCCOMMON"
integer natoms,nangmax,ncontrmax,nprimmax,atnum(*),iatoms,j,icontr
character*2 atsymbol(natoms)
integer nang(natoms),ncontr(0:nangmax,natoms),iang,iprim
integer nprim(0:nangmax,natoms),molperm(*),i,k,l,m,icart,ii,ibasis
real*8 gexp(nprimmax,0:nangmax,natoms),coord(3,natoms)
real*8 gcoef(nprimmax,ncontrmax,0:nangmax,natoms)
character*1 angsym
character*10 cscr
logical cartg
C Write COORD.xyz file
open(moldenfile,file='COORD.xyz')
m=0
do iatoms=1,natoms
if(atnum(iatoms).gt.0) m=m+1
enddo
write(cscr,"(i10)") m
write(moldenfile,"(a10)") adjustl(cscr)
write(moldenfile,*)
do iatoms=1,natoms
if(atnum(iatoms).gt.0) write(moldenfile,"(a2,3f20.10)")
$adjustl(atsymbol(iatoms)),(angtobohr*coord(j,iatoms),j=1,3)
enddo
write(moldenfile,*)
close(moldenfile)
C MOLDEN convention
C 5D: D 0, D+1, D-1, D+2, D-2
C 6D: xx, yy, zz, xy, xz, yz
C 7F: F 0, F+1, F-1, F+2, F-2, F+3, F-3
C 10F: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
C 9G: G 0, G+1, G-1, G+2, G-2, G+3, G-3, G+4, G-4
C 15G: xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
C xxyy xxzz yyzz xxyz yyxz zzxy
C Generate permutation vector to follow the above convention
ibasis=0
ii=0
do iatoms=1,natoms
do iang=0,nang(iatoms)
do icontr=1,ncontr(iang,iatoms)
icart=0
if(cartg) then
do i=0,iang
do k=0,iang-i
m=iang-i-k
icart=icart+1
ibasis=ibasis+1
molperm(ibasis)=ibasis
enddo
enddo
k=ibasis-icart
if(iang.eq.1) then !p functions z,y,x -> x,y,z
molperm(k+1)=k+3
molperm(k+2)=k+2
molperm(k+3)=k+1
else if(iang.eq.2) then
molperm(k+1)=k+6
molperm(k+2)=k+3
molperm(k+3)=k+1
molperm(k+4)=k+5
molperm(k+5)=k+4
molperm(k+6)=k+2
else if(iang.eq.3) then
molperm(k+1)=k+10
molperm(k+2)=k+4
molperm(k+3)=k+1
molperm(k+4)=k+7
molperm(k+5)=k+9
molperm(k+6)=k+8
molperm(k+7)=k+5
molperm(k+8)=k+2
molperm(k+9)=k+3
molperm(k+10)=k+6
else if(iang.eq.4) then
molperm(k+1)=k+15
molperm(k+2)=k+5
molperm(k+3)=k+1
molperm(k+4)=k+14
molperm(k+5)=k+13
molperm(k+6)=k+9
molperm(k+7)=k+4
molperm(k+8)=k+6
molperm(k+9)=k+2
molperm(k+10)=k+12
molperm(k+11)=k+10
molperm(k+12)=k+3
molperm(k+13)=k+11
molperm(k+14)=k+8
molperm(k+15)=k+7
endif
else
do l=-iang,iang
ibasis=ibasis+1
molperm(ibasis)=ibasis
enddo
if(iang.eq.1) then
k=molperm(ibasis)
molperm(ibasis)=molperm(ibasis-2)
molperm(ibasis-2)=k
else if(iang.ge.2) then
k=ibasis-2*iang
molperm(k)=ibasis-iang
do l=1,iang
k=k+1
molperm(k)=ibasis-iang+l
k=k+1
molperm(k)=ibasis-iang-l
enddo
endif
endif
enddo
enddo !iang
enddo !iatoms
C Write MOLDEN file
open(moldenfile,file='MOLDEN.perm')
write(moldenfile,'(100000i7)') (molperm(j),j=1,nbasis)
close(moldenfile)
open(moldenfile,file='MOLDEN')
write(moldenfile,"('[Molden Format]')")
write(moldenfile,"('[ATOMS] AU')")
do iatoms=1,natoms
write(moldenfile,"(a3,2i5,3f20.10)")
$atsymbol(iatoms),iatoms,atnum(iatoms),(coord(j,iatoms),j=1,3)
enddo
write(moldenfile,"('[GTO]')")
do iatoms=1,natoms
write(moldenfile,"(2i4,i2)") iatoms,0
do iang=0,nang(iatoms)
do icontr=1,ncontr(iang,iatoms)
write(moldenfile,"(a2,i5,f5.2)") angsym(iang)
$,nprim(iang,iatoms),1.d0
do iprim=1,nprim(iang,iatoms)
write(moldenfile,"(2d18.10)")
$gexp(iprim,iang,iatoms),gcoef(iprim,icontr,iang,iatoms)
enddo
enddo
enddo
write(moldenfile,*)
enddo
if(.not.cartg) then
write(moldenfile,"('[5D]')")
write(moldenfile,"('[7F]')")
write(moldenfile,"('[9G]')")
endif
write(moldenfile,"('[MO]')")
close(moldenfile)
C
return
end
C
************************************************************************
subroutine embedat(natoms,minpfile,ind,iout,kywd,nkywd)
************************************************************************
* Read atoms of embedded subsystem
************************************************************************
implicit none
integer natoms,m,n,minpfile,iout,i,ii,jj,ind(natoms),nkywd,nn
character*1 line2(16),line3(max(512,6*natoms))
character*8 embed,kywd
character*16 line1
equivalence(line1,line2)
C Read atoms of embedded subsystem
call getkeym(kywd,nkywd,embed,8)
nn=1
if(trim(kywd).eq.'cvs') then
backspace(minpfile)
nn=nkywd+2
endif
n=0
line3=' '
do
n=n+1
read(minpfile,'(a)',advance='no',eor=1234,end=1234) line3(n)
enddo
1234 continue
n=nn
ind=0
do while(line3(n).ne.' ')
m=n
do while(line3(n).ne.','.and.line3(n).ne.'-'
$.and.line3(n).ne.' ')
n=n+1
enddo
line1=' '
do i=m,n-1
line2(i-m+1)=line3(i)
enddo
read(line1,*) ii
jj=ii
if(line3(n).eq.'-') then
n=n+1
m=n
do while(line3(n).ne.','.and.line3(n).ne.' ')
n=n+1
enddo
line1=' '
do i=m,n-1
line2(i-m+1)=line3(i)
enddo
read(line1,*) jj
endif
if(jj.gt.natoms) then
write(iout,*) 'Invalid serial number for '//trim(kywd)//'!'
call mrccend(1)
endif
ind(ii:jj)=1
n=n+1
enddo
C
return
end
C
************************************************************************
subroutine readlinelist(nquant,minpfile,ind,iout,quant)
************************************************************************
* Read quantities list from a line
************************************************************************
implicit none
! intent(in)
integer, intent(in) :: nquant
integer, intent(in) :: minpfile
integer, intent(in) :: iout
integer, intent(out) :: ind(nquant)
! local
character(len=max(512,6*nquant)) :: line
character(len=1) :: line2(16)
character(len=8) :: quant
character(len=16) :: line1
integer :: istat
integer :: i,ii,jj,m,n,kk
integer :: iquant
equivalence(line1,line2)
iquant=0
n=0
line=''
read(minpfile,'(a)',iostat=istat) line
n=1
ind=0
do while(line(n:n).ne.' ')
m=n
do while(line(n:n).ne.','.and.line(n:n).ne.'-'
$.and.line(n:n).ne.' ')
n=n+1
enddo
line1=' '
do i=m,n-1
line2(i-m+1)=line(i:i)
enddo
read(line1,*) ii
if(quant.eq.'mmserial') then
iquant=iquant+1
ind(iquant)=ii
endif
jj=ii
if(line(n:n).eq.'-') then
n=n+1
m=n
do while(line(n:n).ne.','.and.line(n:n).ne.' ')
n=n+1
enddo
line1=' '
do i=m,n-1
line2(i-m+1)=line(i:i)
enddo
read(line1,*) jj
endif
if(jj.gt.nquant.and.quant.ne.'mmserial') then
if(quant.eq.'atoms ') then
write(iout,*) 'Invalid serial number for atoms!'
else if(quant.eq.'mos ') then
write(iout,*) 'Invalid serial number for MOs!'
endif
call mrccend(1)
endif
if(quant.ne.'mmserial') then
ind(ii:jj)=1
else
do kk=ii+1,jj
iquant=iquant+1
ind(iquant)=kk
enddo
endif
n=n+1
enddo
if(istat.ne.0) then
write(iout,'(a)') 'Cannot read a list of '//trim(adjustl(quant))
$//' from MINP file.'
write(iout,'(a)')
$' The error is detected @ readlinelist (combin.f)'
call mrccend(1)
endif
C
return
end
C
************************************************************************
logical function quadratic(a,b,c,r1,r2)
************************************************************************
* Calculate the roots of a general quadratic equation
************************************************************************
implicit none
real*8 a,b,c,r1,r2,disc
C
disc=b*b-4.d0*a*c
if(disc.lt.0.d0) then
quadratic=.false.
else
quadratic=.true.
disc=dsqrt(disc)
r1=(-b+disc)/(2.d0*a)
r2=(-b-disc)/(2.d0*a)
endif
C
return
end
C
************************************************************************
subroutine prtenergc(meth,ene,erefc,locnoc)
************************************************************************
* Print energies *
************************************************************************
#include "MRCCCOMMON"
real*8 ene,erefc
integer locnoc
character*15 beg
character*(*) meth
C
call timer
#if defined(MPI)
if(.not.master_thread .and. locnoc.eq.0) return ! need to print it if lccoporder=lccfirst & MPI
#endif
c if(log1) then
c beg=' Total ' // meth
c end=' energy [au]:'
c call trimloc(beg,end)
c write(iout,"(a31,f22.12)") beg,ene
c endif
beg=' '
beg=meth
if(locnoc.gt.0) then
write(ifcfile,7596)
$'ENERGY ','REF ',1,1,1,erefc,0.d0,0.d0
write(ifcfile,7596)
$'ENERGY ',beg,1,1,1,ene,tim-timold,wct-wctold
else if(trim(beg).eq.'CCSD(T)') then
write(ifcfile,7596)
$'ENERGY ',beg,1,1,1,ene+erefc,tim-timold,wct-wctold
else
write(ifcfile,7596)
$'ENERGY ',beg,1,1,1,ene,tim-timold,wct-wctold
endif
timold=tim
wctold=wct
7596 format(a8,1x,a15,1x,3(i2,1x),1pe23.15,2(1pe15.8))
C
return
end
C
************************************************************************
subroutine genprojmx(dens,cmo,nb,nmo,ss,tol,scmo,smo,redunant,scr,
$mode)
************************************************************************
* construct the mx projecting onto a possibly overlapping cmo basis
* if redundant, linear dependency of the cmo basis is removed via discarding
* S_cmo eigenvectors with lower than tol scmo eigenvalues
* mode=1 make full projector mx
* mode=2 make only the invrese of the cmo overlap mx
* mode=3 make the inverse of the input overlap mx (provided in smo)
************************************************************************
implicit none
integer nb,nmo,mode
real*8 dens(nb,nb),cmo(nb,nmo),ss(nb,nb),tol
real*8 scmo(nb,nmo),smo(nmo,nmo),scr(*)
logical redunant
c internal
integer iisyev,nred,j,i
integer*4 isyev
equivalence(isyev,iisyev) !For Intel
C
if (mode.eq.1.or.mode.eq.2) then
c overlap mx of the cmo basis goes to Smo
call dsymm('l','l',nb,nmo,1.d0,ss,nb,cmo,nb,0.d0,scmo,nb)
call dgemm('t','n',nmo,nmo,nb,1.d0,cmo,nb,scmo,nb,0.d0,
$smo,nmo)
endif
c invert Smo
if (redunant) then
call dsyev('V','U',nmo,smo,nmo,scr,scr(1+nmo),20*nmo,isyev)
if(isyev.ne.0) then
write(*,*) 'Error at inversion in genprojmx'
call mrccend(1)
endif
nred=0
c write(*,*) 'S_paop eigvals'
c write(*,'(10000es14.6)') scr(1:nmo)
do while(scr(nred+1).lt.tol)
nred=nred+1
c write(*,*) 'drop redundant cmo',nred,scr(nred)
if (scr(nred+1).lt.-tol) then
write(*,*) 'warning <0 S_mo eigval',nred,scr(nred)
endif
enddo
c if (nred.gt.0)
c $ write(*,*) 'drop',nred,'redundant cmo in genprojmx',scr(nred)
do j=nred+1,nmo
call dscal(nmo,1.d0/dsqrt(scr(j)),smo(1,j),1)
enddo
call dsyrk('u','n',nmo,nmo-nred,1.d0,smo(1,nred+1),nmo,0.d0,
$scr,nmo)
call dcopy(nmo**2,scr,1,smo,1)
call filllo(smo,nmo)
else
call syminvpd(smo,nmo,nmo)
call fillup(smo,nmo)
endif
if (mode.eq.1) then
c return the mx of the projector: cmo*Smo^-1*cmo^T
call dsymm('r','u',nb,nmo,1.d0,smo,nmo,cmo,nb,0.d0,scmo,nb)
call dgemm('n','t',nb,nb,nmo,1.d0,scmo,nb,cmo,nb,0.d0,dens,nb)
endif
c
return
end
c
************************************************************************
subroutine normalize(nb,nmo,cmo,scmo,ss,havescmo,norm,
& normalize_scmo)
************************************************************************
* Normalize the input orbitals in cmo, scmo=overlap*cmo, ss: overlap mx of the basis functions
************************************************************************
implicit none
integer nb,nmo,i
real*8 cmo(nb,nmo),ss(nb,nb)
real*8 scmo(nb,nmo),tmp,ddot,norm
logical havescmo,normalize_scmo
c
if (.not.havescmo) then
call dsymm('l','l',nb,nmo,1.d0,ss,nb,cmo,nb,0.d0,scmo,nb)
endif
c
do i=1,nmo
tmp=ddot(nb,cmo(1,i),1,scmo(1,i),1)
tmp=dsqrt(norm/tmp)
call dscal(nb,tmp,cmo(1,i),1)
if (normalize_scmo) call dscal(nb,tmp,scmo(1,i),1)
enddo
C
return
end
c
subroutine canonicalize(nb,nmo,fmat,cmo,pscan,scr,ev)
************************************************************************
* Create pseudo canonical orbitals of the subspace spanned by the input orbitals in cmo
* fmat: Fock mx
************************************************************************
implicit none
integer nb,nmo,i,iisyev,j
real*8 cmo(nb,nmo),fmat(nb,nb),ev(nmo),pscan(nmo,nmo),scr(*)
integer*4 isyev
equivalence(isyev,iisyev) !For Intel
c
c Fock mx of the cmo basis
call dsymm('l','l',nb,nmo,1.d0,fmat,nb,cmo,nb,0.d0,scr,nb)
call dgemm('t','n',nmo,nmo,nb,1.d0,cmo,nb,scr,nb,0.d0,pscan,nmo)
c
call dsyev('V','U',nmo,pscan,nmo,ev,scr,20*nmo,isyev)
if(isyev.ne.0) then
write(*,*) 'Error in subroutine canonicalize'
call mrccend(1)
endif
c
return
end
c
subroutine lwdncanOG(nb,nmo,ss,cmo,scmo,scr,tol,havescmo,lindep,ev
$,smo,ndrop,overwrite,havesmo,ogcmo)
************************************************************************
* Orthogonalize and remove linearly dependent cmos, if needed
* scmo=overlap*cmo, ss: overlap mx of the basis functions
* if havesmo then the overlap is not computed
* if havescmo then AO overlap x input mo coeff is not computed
************************************************************************
implicit none
integer nb,nmo,i,iisyev,j,ndrop
real*8 cmo(nb,nmo),ss(nb,nb),scmo(nb,*),scr(*),ev(nmo),smo(nmo,*)
real*8 tol,ogcmo(nb,nmo)
integer*4 isyev
equivalence(isyev,iisyev) !For Intel
logical havescmo,lindep,overwrite,havesmo
c
if (.not.havesmo) then
if (.not.havescmo) then
call dsymm('l','l',nb,nmo,1.d0,ss,nb,cmo,nb,0.d0,scmo,nb)
endif
call dgemm('t','n',nmo,nmo,nb,1.d0,cmo,nb,scmo,nb,0.d0,smo,nmo)
endif
c write(*,"('norms in lwdncanOG:',10000es14.6)")(smo(j,j),j=1,nmo)
scr(1)=1.d0
do j=1,nmo
if (smo(j,j).lt.scr(1)) scr(1)=smo(j,j)
if (smo(j,j).lt.1d-4) write(*,*) 'small norm',j,smo(j,j)
enddo
c write(*,"('min norm in lwdncanOG:',10es14.6)")scr(1)
C Remove linearly dependent orbitals and orthogonalize
call dsyev('V','U',nmo,smo,nmo,ev,scr,20*nmo,isyev)
if(isyev.ne.0) then
write(*,*) 'Error in subroutine lwdncanOG'
call mrccend(1)
endif
ndrop=0
if (lindep) then
do while(ev(ndrop+1).lt.tol)
ndrop=ndrop+1
c write(*,*) 'drop redundant orb', ndrop,ev(ndrop)
enddo
endif
c write(*,"('next S ev in lwdncanOG:',1000es14.6)")
c $ev(ndrop+1:min(ndrop+10,nmo))
do j=ndrop+1,nmo
call dscal(nmo,1.d0/dsqrt(ev(j)),smo(1,j),1)
enddo
call dgemm('n','n',nb,nmo-ndrop,nmo,1.d0,cmo,nb,smo(1,ndrop+1),nmo
$,0.d0,ogcmo,nb)
if (overwrite) call dcopy(nb*(nmo-ndrop),ogcmo,1,cmo,1)
c
return
end
c
c {{{ subroutine dfillzeroup
subroutine dfillzeroup(a,n)
************************************************************************
* Fill up the upper triangle of a matrix with zeros
************************************************************************
implicit none
integer n,i,j
real*8 a(n,n)
C
do i=1,n
do j=i+1,n
a(i,j)=0.d0
enddo
enddo
C
return
end
c }}}
c
c {{{ SUBROUTINE GENTRANS
SUBROUTINE GENTRANS(A,RD,M,N,MN,B,MOVE,IWRK,IOK)
C IN-PLACE TRANSFORM OF A(N ARRAY OF) RECTANGULAR OR SQUARE MATRIX(S)
C INPUT:
C A IS A MATRIX WITH RD ROW DIMENSION AND MN=M*N COLUMN DIMENSION
C B IS A BUFFER FOR THE ELEMENT INTERCHANGE WITH RD DIMENSION
C THE M BY N MATRIX OF A(RD,.) IS TO BE TRANSPOSED FOR EACH RD VALUE
C MOVE IS A ONE-DIMENSIONAL ARRAY OF LENGTH IWRK
C VALUE IWRK=(M+N)/2 IS RECOMMENDED.
C IOK INDICATES THE SUCCESS OR FAILURE OF THE ROUTINE.
C NORMAL RETURN IOK=0
C ERRORS IOK=-1, MN NOT EQUAL TO M*N.
C IOK=-2, IWRK NEGATIVE OR ZERO.
C IOK.GT.0, (SHOULD NEVER OCCUR). IN THIS CASE
C WE SET IOK EQUAL TO THE FINAL VALUE OF I WHEN THE SEARCH
C IS COMPLETED BUT SOME LOOPS HAVE NOT BEEN MOVED.
C
C BASED ON THE ALGORITH IN COMM. ACM, VOL. 13, NO. 05, P. 324. AND NETLIB.ORG/TOMS/380
IMPLICIT NONE
INTEGER IWRK,MOVE(IWRK),RD,MN,M,N,NCOUNT,M2,I,IB,IA,K,KMI,MAX,I1
INTEGER I2,IOK,N1,J,J1
real*8 A(RD,MN),B(RD)
C CHECK ARGUMENTS AND INITIALISE
IF(M.LT.2.OR.N.LT.2) GO TO 60
IF(MN.NE.M*N) GO TO 92
IF(IWRK.LT.1) GO TO 93
IF(M.EQ.N) GO TO 70
NCOUNT=2
M2=M-2
DO 10 I=1,IWRK
10 MOVE(I)=0
IF(M2.LT.1) GO TO 12
C COUNT NUMBER,NCOUNT, OF SINGLE POINTS.
DO 11 IA=1,M2
IB=IA*(N-1)/(M-1)
IF(IA*(N-1).NE.IB*(M-1)) GO TO 11
NCOUNT=NCOUNT+1
I=IA*N+IB
IF(I.GT.IWRK) GO TO 11
MOVE(I)=1
11 CONTINUE
C SET INITIAL VALUES FOR SEARCH.
12 K=MN-1
KMI=K-1
MAX=MN
I=1
C AT LEAST ONE LOOP MUST BE RE-ARRANGED.
GO TO 30
C SEARCH FOR LOOPS TO BE REARRANGED.
20 MAX=K-I
I=I+1
KMI=K-I
IF(I.GT.MAX) GO TO 90
IF(I.GT.IWRK) GO TO 21
IF(MOVE(I).LT.1) GO TO 30
GO TO 20
21 IF(I.EQ.M*I-K*(I/N)) GO TO 20
I1=I
22 I2=M*I1-K*(I1/N)
IF(I2.LE.I .OR. I2.GE.MAX) GO TO 23
I1=I2
GO TO 22
23 IF(I2.NE.I) GO TO 20
C REARRANGE ELEMENTS OF A LOOP.
30 I1=I
31 B(1:RD)=A(1:RD,I1+1)
32 I2=M*I1-K*(I1/N)
IF(I1.LE.IWRK) MOVE(I1)=2
33 NCOUNT=NCOUNT+1
IF(I2.EQ.I .OR. I2.GE.KMI) GO TO 35
34 A(1:RD,I1+1)=A(1:RD,I2+1)
I1=I2
GO TO 32
35 IF(MAX.EQ.KMI .OR. I2.EQ.I) GO TO 41
MAX=KMI
GO TO 34
C TEST FOR SYMMETRIC PAIR OF LOOPS.
41 A(1:RD,I1+1)=B(1:RD)
IF(NCOUNT.GE.MN) GO TO 60
IF(I2.EQ.MAX .OR. MAX.EQ.KMI) GO TO 20
MAX=KMI
I1=MAX
GO TO 31
C NORMAL RETURN.
60 IOK=0
RETURN
C IF MATRIX IS SQUARE, EXCHANGE ELEMENTS A(I,J) AND A(J,I).
70 N1=N-1
DO 71 I=1,N1
J1=I+1
DO 71 J=J1,N
I1=I+(J-1)*N
I2=J+(I-1)*M
B(1:RD)=A(1:RD,I1)
A(1:RD,I1)=A(1:RD,I2)
A(1:RD,I2)=B(1:RD)
71 CONTINUE
GO TO 60
C ERROR RETURNS.
90 IOK=I
91 RETURN
92 IOK=-1
GO TO 91
93 IOK=-2
GO TO 91
END
c }}}
c
c {{{ subroutine relocatemx
subroutine relocatemx(m1,m2,drow,dcol,backward)
************************************************************************
* copy the mx in m1 to m2 in memory elementwise
* hence m1 and m2 may be overlapping in memory
************************************************************************
#if defined (Intel)
c !DIR$ NOOPTIMIZE
!DIR$ OPTIMIZE:1
#endif
implicit none
integer i,j,drow,dcol
real*8 m1(drow,dcol),m2(drow,dcol)
logical backward
C
if (backward) then
c loc(m1) > loc(m2)
!$OMP SINGLE
do i=1,dcol
do j=1,drow
m2(j,i)=m1(j,i)
enddo
enddo
!$OMP END SINGLE
else
c loc(m1) < loc(m2)
!$OMP SINGLE
do i=dcol,1,-1
do j=drow,1,-1
m2(j,i)=m1(j,i)
enddo
enddo
!$OMP END SINGLE
endif
C
return
end
c }}}
c
c {{{ function localcc15p
logical function localcc15p(localcc)
************************************************************************
* T: if localcc >=2015 , F: otherwise
************************************************************************
implicit none
character*4 localcc
c
localcc15p=localcc.eq.'2015'.or.localcc.eq.'2016'.or.
$ localcc.eq.'2018'.or.localcc.eq.'2021'.or.localcc.eq.'2022'
C
return
end
c }}}
c
c {{{ function localccXp
logical function localccXp(localcc,year)
************************************************************************
* T: if localcc >=year , F: otherwise
************************************************************************
implicit none
character*4 localcc
integer year,localcc2int
c
if (localcc.eq.'off ') then
localccXp=.false.
elseif (localcc.eq.'on '.or.localcc.eq.' ') then
write(*,*)'Error: unsupported localcc value in localccXp',localcc
call mrccend(1)
else
read(localcc,*) localcc2int
localccXp=localcc2int.ge.year
endif
C
return
end
c }}}
c
c {{{ lstartwl and lstartwlno
************************************************************************
logical function lstartwl(calc)
************************************************************************
* Check if calc starts with l, true if a local correlation computation is requested
************************************************************************
implicit none
character(len=32) calc,xcalc
character*1 calc1(32)
equivalence(xcalc,calc1)
c
xcalc=calc
lstartwl=calc1(1).eq.'l'.and.(calc1(2).ne.'d'.or.calc1(3).ne.'a')
$ .and.(calc1(2).ne.'c'.or.calc1(3).ne.'-')
return
end
C
************************************************************************
logical function lstartwlno(calc)
************************************************************************
* Check if calc starts with LNO-, true if a local correlation computation is requested
************************************************************************
implicit none
character(len=32) calc,xcalc
character*1 calc1(32)
equivalence(xcalc,calc1)
c
xcalc=calc
lstartwlno=calc1(1).eq.'l'.and.calc1(2).eq.'n'.and.calc1(3).eq.'o'
$.and.calc1(4).eq.'-'
return
end
C }}}
c
************************************************************************
subroutine removeprefix(calc1,clen,nchartoremove)
************************************************************************
* remove the first nchartoremove character of calc, e.g. for DF-, RI-, LNO-, L prefixes
************************************************************************
implicit none
integer clen,nchartoremove,i
character*1 calc1(clen)
c
do i=1,clen-nchartoremove
calc1(i)=calc1(i+nchartoremove)
enddo
do i=clen-nchartoremove+1,clen
calc1(i)=' '
enddo
c
return
end
************************************************************************
subroutine basislocation(cdmrcc,cdmrcc1)
************************************************************************
* determine the folder containing the basis set information
************************************************************************
#include "MRCCCOMMON"
integer*4 istat
character*512 cdmrcc
character*1 cdmrcc1(512)
character*128 command
integer i
#if !defined (gfortran) && !defined (G95)
integer*4 system
external system
#endif
c
cdmrcc1=' '
command=trim("which dmrcc > mrccjunk1") // char(0)
istat=system(command)
if(istat.eq.0) then
open(unit=scrfile1,file='mrccjunk1')
read(scrfile1,'(a512)') cdmrcc
cdmrcc=adjustl(cdmrcc)
close(scrfile1,status='delete')
i=507
do while(cdmrcc1(i ).ne.'d'.or.
$ cdmrcc1(i+1).ne.'m'.or.
$ cdmrcc1(i+2).ne.'r'.or.
$ cdmrcc1(i+3).ne.'c'.or.
$ cdmrcc1(i+4).ne.'c')
i=i-1
enddo
cdmrcc1(i )='B'
cdmrcc1(i+1)='A'
cdmrcc1(i+2)='S'
cdmrcc1(i+3)='I'
cdmrcc1(i+4)='S'
cdmrcc1(i+5)='/'
cdmrcc1(i+6:512)=' '
else
write(*,*) 'Error in determining PATH of BASIS folder',istat
write(*,*) ' Warning: assuming ./BASIS'
cdmrcc1(1)='.'
cdmrcc1(2)='/'
i=3
cdmrcc1(i )='B'
cdmrcc1(i+1)='A'
cdmrcc1(i+2)='S'
cdmrcc1(i+3)='I'
cdmrcc1(i+4)='S'
cdmrcc1(i+5)='/'
cdmrcc1(i+6:512)=' '
endif
c
return
end
c
************************************************************************
function iname(ch)
implicit none
integer iname
character*4 ch
if(ch.eq."abcd")iname=0
if(ch.eq."ijkl")iname=1
if(ch.eq."abci")iname=2
if(ch.eq."aijk")iname=3
if(ch.eq."abij")iname=4
if(ch.eq."aibj")iname=5
if(ch.eq."iajk")iname=6
if(ch.eq."abic")iname=7
if(ch.eq."aibc")iname=8
if(ch.eq."iabc")iname=9
if(ch.eq."ijak")iname=10
if(ch.eq."ijka")iname=11
if(ch.eq."iajb")iname=12
if(ch.eq."aijb")iname=13
if(ch.eq."ijab")iname=14
if(ch.eq."iabj")iname=15
end
c
************************************************************************
subroutine write56(nbmax,iface,ccalc,domain_is_cs,lis_nal,lis_nbe)
************************************************************************
* Write fort.56 for ovirt and mrcc
************************************************************************
#include "MRCCCOMMON"
c taken from MRCCCOMMON: ncore,op,calc,nocc,rest,dens,conver,isym,diag,
c iclsh,ihf,ndoub,nsing,ntrip,nacto,nactv,maxex,ptfreq,sacc,dboc,
c maxmem,maxcor,idnint
c minpfile,eifile,iout
integer i,j,k,m,n,ii,jj,iitol,nbmax,occvec(nbmax),iacto,iactv,nco
integer iocc,lis_nal,lis_nbe,nact,ndoc
real*8 eps,rr
character*1 xcalc1(16),chs,line2(16),line3(512)
character*2 ch2
character*4 cmult,cctol,symm,localcc,nstate,nnsing,nntrip,ovirt
character*4 ch4
character*5 scftype
character*6 cdiag,core
character*8 active,refdet,iface,embed
character*13 rohftype
character*16 ccalc,xcalc,mem,orblocc,orbloco,orblocv,orbloce,line1
character*512 line
logical domain_is_cs
equivalence(line,line3)
equivalence(line1,line2)
equivalence(xcalc,xcalc1)
C
call getkey('embed',5,embed,8)
open(minpfile,file='MINP')
call getkey('core',4,core,6)
nco=ncore
if(core.eq.'corr ') nco=0
if(core.eq.'corr '.and.embed.eq.'huzinaga') nco=nfroz
C ex.lev. & calc
c call getkey('calc',4,ccalc,16)
op=0
calc=1
if(ccalc.eq.'fci ') then
op=nocc-2*nco
calc=0
endif
xcalc=ccalc
if(xcalc1(1).eq.'c'.and.xcalc1(2).eq.'i') calc=0
xcalc1(1)=' '
xcalc1(2)=' '
do i=3,15
if(xcalc1(i).eq.'/') calc=9
if(xcalc1(i).eq.'_') calc=4
if(xcalc1(i).eq.'-') then
if(xcalc1(i+1).eq.'3') calc=8
if(xcalc1(i+1).eq.'1') then
if(xcalc1(i+2).eq.'b') then
calc=6
else
calc=5
endif
endif
do j=i,16
xcalc1(j)=' '
enddo
endif
enddo
if(xcalc.eq.' s ') op=1
if(xcalc.eq.' sd ') op=2
if(xcalc.eq.' sdt ') op=3
if(xcalc.eq.' sdtq ') op=4
if(xcalc.eq.' sdtqp ') op=5
jj=0
do i=3,16
ii=ichar(xcalc1(i))
if((ii.lt.49.or.ii.gt.57).and.xcalc1(i).ne.' ') jj=jj+1
enddo
if(jj.eq.0) then
calc=7
read(xcalc,*) op
endif
if(op.eq.0) then
if(calc.le.1) then
do i=3,15
if(xcalc1(i).eq.')'.and.xcalc1(i+1).eq.'(') calc=3
if(xcalc1(i).eq.')'.and.xcalc1(i+1).eq.'[') calc=2
enddo
endif
i=16
do while(xcalc1(i).ne.')'.and.xcalc1(i).ne.']')
i=i-1
enddo
j=i-2
chs='('
if(xcalc1(i).eq.']') chs='['
do while(xcalc1(j).ne.chs)
j=j-1
enddo
do ii=1,j
xcalc1(ii)=' '
enddo
do ii=i,16
xcalc1(ii)=' '
enddo
if(xcalc1(j+1).eq.'d') then
op=2
else if(xcalc1(j+1).eq.'t') then
op=3
if(calc.le.1) calc=3
else if(xcalc1(j+1).eq.'q') then
op=4
if(calc.le.1) calc=3
else if(xcalc1(j+1).eq.'p') then
op=5
if(calc.le.1) calc=3
else
read(xcalc,*) op
endif
if(calc.eq.3.and.chs.eq.'[') calc=2
endif
if(op.gt.nal+nbe-2*nco) then
op=nal+nbe-2*nco
if(calc.gt.1) calc=1
endif
C rest & fort.57
call getkey('rest',4,ch4,4)
read(ch4,*) rest
if(rest.ge.3) then
call getkeym('rest',4,ch4,4)
read(minpfile,'(512a1)') line3
open(eifile,file='fort.57',status='unknown')
rewind(eifile)
write(eifile,*) trim(line)
close(eifile)
else
open(eifile,file='fort.57',status='unknown')
close(eifile,status='delete')
endif
call getkey('ciguess',7,ch4,4)
if(ch4.eq.'on ') then
call getkeym('ciguess',7,ch4,4)
c call ishell("sed '1,/ciguess/d' MINP |grep -iv [a-z] > fort.57")
open(eifile,file='fort.57',status='unknown')
do
read(minpfile,'(512a1)',end=8599) line3
if(ichar(line3(1)).lt.48.or.ichar(line3(1)).gt.57) exit
write(eifile,*) trim(line)
enddo
8599 continue
close(eifile)
endif
C dens
call getkey('dens',4,ch4,4)
read(ch4,*) dens
if(mod(dens,2).eq.1) dens=dens+1
if(iface.eq.'columbus'.or.iface.eq.'none ') dens=-dens
C conver
conver=0
C symm
call getkey('symm',4,symm,4)
call getkey('localcc',7,localcc,4)
call getkey('orblocc',7,orblocc,16)
call getkey('orbloco',7,orbloco,16)
call getkey('orblocv',7,orblocv,16)
call getkey('orbloce',7,orbloce,16)
if(trim(orblocc).ne.'off'.or.
$ trim(orbloco).ne.'off'.or.
$ trim(orblocv).ne.'off'.or.
$ trim(orbloce).ne.'off') then
isym=0
else if(symm.eq.' '.or.localcc.ne.'off ') then
isym=1
else
read(symm,*) isym
endif
C diag
call getkey('diag',4,cdiag,6)
if(cdiag.eq.'olsen ') then
diag=1
else if(cdiag.eq.'follow') then
diag=2
else
diag=0
endif
C CS
call getkey('mult',4,cmult,4)
if(cmult.eq.'1 ') then
iclsh=1
else
iclsh=0
endif
if(iface.eq.'dirac ') iclsh=0
C spatial
call getkey('scftype',7,scftype,5)
call getkey('rohftype',8,rohftype,13)
if((scftype.eq.'uhf '.or.(scftype.eq.'rohf '.and.
$rohftype.eq.'semicanonical')).and..not.domain_is_cs) then
k=0
else
k=1
endif
if(iface.eq.'dirac ') k=1
C HF
if(scftype.eq.'rhf '.or.scftype.eq.'uhf ') then
ihf=1
else
ihf=0
endif
if(trim(orblocc).ne.'off'.or.
$ trim(orbloco).ne.'off'.or.
$ trim(orblocv).ne.'off'.or.
c $ trim(orbloce).ne.'off'.or.
$ embed.eq.'scl '.or.embed.eq.'fdm ') ihf=0
if(localcc.ne.'off '.and.embed.ne.'scl '.and.
$ embed.ne.'fdm ') ihf=1
C nsing, ntrip, & ndoub
call getkey('nstate',6,nstate,4)
read(nstate,*) ndoub
call getkey('nsing' ,5,nnsing,4)
read(nnsing,*) nsing
call getkey('ntrip' ,5,nntrip,4)
read(nntrip,*) ntrip
if(iclsh.eq.1) then !.and.k.eq.1) then
if(nsing.eq.0) nsing=ndoub
if(calc.ge.1) nsing=max(1,nsing)
ndoub=0
else
ndoub=max(ndoub,nsing+ntrip)
nsing=0
ntrip=0
endif
C nacto
call getkey('nacto',5,ch4,4)
read(ch4,*) nacto
C nactv
call getkey('nactv',5,ch4,4)
read(ch4,*) nactv
C tol
call getkey('cctol',5,cctol,4)
read(cctol,*) iitol
C maxex
call getkey('maxex',5,ch4,4)
read(ch4,*) maxex
!!! !HB
if((maxex.eq.0.and.embed.eq.'huzinaga').or.nvfroz.ne.0) then
maxex=op
endif
!!! !HB
C sacc
sacc=0
C ptfreq
call getkey('ptfreq',6,line1,16)
read(line1,*) ptfreq
C dboc
dboc=0
call getkey('dboc',4,ch4,4)
if(ch4.eq.'on ') dboc=1
C mem
call getkey('mem',3,line1,16)
if(line1.ne.' ') then
i=1
do while(line2(i).ne.' '.and.i.lt.16)
i=i+1
enddo
call lowercase(line2(i-2),ch2,2)
if(ch2.ne.'mb'.and.ch2.ne.'gb') then
write(iout,*) 'Unknown memory unit!'
call exit(1)
endif
line2(i-2)=' '
line2(i-1)=' '
read(line1,*) maxmem
else
ch2='mb'
maxmem=256.d0
endif
if(ch2.eq.'gb') maxmem=1024.d0*maxmem
maxcor=idnint(maxmem)
C locno & eps
locno=0
eps=0.d0
call getkey('eps',3,line1,16)
call getkey('ovirt',5,ovirt,4)
if(ovirt.eq.'mp2 ') locno=1
if(ovirt.eq.'ovos') locno=2
if(ovirt.ne.'off ') then
if(line1 .ne. ' ') then
read(line1,*) eps
else
call getkey('lnoepsv',7,line1,16)
read(line1,*) eps
end if
end if
call getkey('localcc',7,localcc,4)
if(localcc.ne.'off ') then
locno=3
call getkey('lnoepsv',7,line1,16)
read(line1,*) eps
call getkey('lnoepso',7,line1,16)
read(line1,*) rr
eps=min(eps,rr)
endif
C Occupation vector
call getkey('refdet',6,refdet,8)
if(refdet.ne.'none ') then
call getkeym('refdet',6,refdet,8)
if(refdet.eq.'serialno') then
call ifillzero(occvec,nbasis-nco)
do iocc=2,-1,-1
if(iocc.ne.0) then
read(minpfile,'(512a1)') line3
line=adjustl(line)
n=1
do while(line3(n).ne.' ')
m=n
do while(line3(n).ne.','.and.line3(n).ne.'-'
$.and.line3(n).ne.' ')
n=n+1
enddo
line1=' '
do i=m,n-1
line2(i-m+1)=line3(i)
enddo
read(line1,*) ii
jj=ii
if(line3(n).eq.'-') then
n=n+1
m=n
do while(line3(n).ne.','.and.line3(n).ne.' ')
n=n+1
enddo
line1=' '
do i=m,n-1
line2(i-m+1)=line3(i)
enddo
read(line1,*) jj
endif
if(jj.gt.nbasis-nco) then
write(iout,*)
$'Invalid serial number for keyword refdet!'
call dmrccend(1)
endif
do i=ii,jj
occvec(i)=iocc
enddo
n=n+1
enddo
endif
enddo
else
read(minpfile,*) (occvec(i),i=1,nbasis-nco)
endif
do i=1,nbasis-nco
if(occvec(i).ne.2.or.occvec(i).ne.0) then
iclsh=0
ihf=0
ndoub=ndoub+nsing+ntrip
nsing=0
ntrip=0
exit
endif
enddo
endif
C Options (first line)
open(inpfile,status='unknown',file='fort.56')
rewind(inpfile)
write(inpfile,"(i4,17i7,f9.5,i7,i9,i5,1pe10.3)")
c write(inpfile,"(i4,17i7,f9.5,3i7,1pe10.3)")
c write(inpfile,"(i4,17i6,f9.5,3i5,1pe10.3)")
$op,nsing,ntrip,rest,calc,dens,conver,isym,diag,iclsh,k,ihf,ndoub,
$nacto,nactv,iitol,maxex,sacc,ptfreq,dboc,maxcor,locno,eps
C Keword names (second line)
write(inpfile,"(a138)") 'ex.lev,nsing,ntrip, rest,CC/CI,dens,' //
$'conver,symm, diag, CS ,spatial, HF ,ndoub,nacto,nactv, ' //
$'tol ,maxex, sacc, freq, dboc, mem, locno, eps'
C Occupation vector (third line)
if (localcc.ne.'off '.and.trim(scftype).ne.'rhf') then
write(inpfile,"(1000000i2)") (2,i=1,lis_nbe),
$ (1,i=1,lis_nal-lis_nbe),(0,i=1,nbasis-lis_nal)
endif !localcc.ne.'off'
if(localcc.eq.'off ') then
if(refdet.ne.'none ') then
write(inpfile,"(100000i3)") (occvec(i),i=1,nbasis-nco)
else
if(iface.eq.'none ') then
write(inpfile,"(1000000i2)") (2,i=1,nbe-nco),(1,i=1,nal-nbe),
$(0,i=1,nbasis-nal)
else
write(inpfile,"(1000000i3)") (2,i=1,nbe),(1,i=1,nal-nbe),
$(0,i=1,nbasis-nal)
endif
endif
C Active/inactive oribitals (fourth line)
rewind(inpfile)
read(inpfile,*)
read(inpfile,*)
read(inpfile,*)
if((embed.eq.'huzinaga'.and.nfroz.ne.0).or.nvfroz.ne.0) then
call ifillzero(occvec,nbasis-nco)
occvec(1:nal-nco)=1 ! occupied mos
occvec(nal-nco+nvfroz+1:nbasis-nco)=1 ! virtual mos
write(inpfile,"(100000i2)") (occvec(i),i=1,nbasis-nco)
else
call getkey('active',6,active,8)
if(active.ne.'none ') then
call getkeym('active',6,active,8)
if(active.eq.'serialno') then
call ifillzero(occvec,nbasis-nco)
read(minpfile,'(512a1)') line3
line=adjustl(line)
n=1
do while(line3(n).ne.' ')
m=n
do while(line3(n).ne.','.and.line3(n).ne.'-'
$.and.line3(n).ne.' ')
n=n+1
enddo
line1=' '
do i=m,n-1
line2(i-m+1)=line3(i)
enddo
read(line1,*) ii
jj=ii
if(line3(n).eq.'-') then
n=n+1
m=n
do while(line3(n).ne.','.and.line3(n).ne.' ')
n=n+1
enddo
line1=' '
do i=m,n-1
line2(i-m+1)=line3(i)
enddo
read(line1,*) jj
endif
if(jj.gt.nbasis-nco) then
write(iout,*)'Invalid serial number for active orbitals!'
call dmrccend(1)
endif
do i=ii,jj
occvec(i)=1
enddo
n=n+1
enddo
else
read(minpfile,*) (occvec(i),i=1,nbasis-nco)
endif
write(inpfile,"(100000i2)") (occvec(i),i=1,nbasis-nco)
else if(trim(scftype).eq.'mcscf') then
call nactorb('docc',ndoc)
call nactorb('mact',nact)
write(inpfile,"(100000i2)") (0,i=1,ndoc-nco),(1,i=1,nact),
$(0,i=nact+1,nbasis-nco)
endif
endif
C Maximum number of inactive orbitals (fifth line)
if((embed.eq.'huzinaga'.and.nfroz.ne.0).or.nvfroz.ne.0) then
call ifillzero(occvec,nbasis-nco)
write(inpfile,"(100000i3)") (occvec(i),i=1,op)
else
call getkey('maxact',6,ch4,4)
if(ch4.eq.'on ') then
if(active.eq.'none ') then
backspace(inpfile)
read(inpfile,*) (occvec(i),i=1,nbasis-nco)
iacto=nacto
iactv=nactv
do i=1,nbasis-nco
occvec(i)=iabs(occvec(i))
enddo
ii=nbasis-nco
do while(occvec(ii).eq.0)
ii=ii-1
enddo
jj=ii
do while(occvec(ii).eq.1)
ii=ii-1
iacto=iacto-1
iactv=iactv-1
enddo
do i=jj+1,jj+iactv/2
occvec(i)=1
enddo
do i=ii,ii-iacto/2+1,-1
occvec(i)=1
enddo
do i=1,nbasis-nco
occvec(i)=mod(occvec(i),2)
enddo
write(inpfile,"(100000i2)") (occvec(i),i=1,nbasis-nco)
endif
call getkeym('maxact',6,ch4,4)
read(minpfile,*) (occvec(i),i=1,op)
write(inpfile,"(100000i3)") (occvec(i),i=1,op)
endif
endif
endif ! localcc.eq.off
C
close(minpfile)
close(inpfile)
C
return
end
c
************************************************************************
subroutine sym2aobasistrf(symmat,aomat,n,nir,nfunc,offset,
$sqroffset,scr,scr2,scrfile1)
************************************************************************
* transform symmat from symmerty basis to the original AO basis of dimension n
* output: aomat
************************************************************************
implicit none
integer n,nir,nfunc(nir),sqroffset(nir),offset(nir),scrfile1,i
real*8 scr(*),scr2(*),symmat(*),aomat(*)
c
if(nir.gt.1) then
open(scrfile1,file='SYMTRA',form='unformatted')
call rtdmx(scr2,scr2,scr,scrfile1,n,n)
close(scrfile1)
do i=1,nir
if(nfunc(i).gt.0) then
call dsymm('r','u',n,nfunc(i),1.d0,symmat(sqroffset(i)),
$nfunc(i),scr((offset(i)-1)*n+1),n,0.d0,scr2((offset(i)-1)*n+1),n)
endif
enddo
call dgemm('n','t',n,n,n,1.d0,scr2,n,scr,n,0.d0,aomat,n)
endif
c
return
end
************************************************************************
subroutine qsi(iwork,n,a)
************************************************************************
* Integer quick sort
************************************************************************
implicit none
integer iwork(*),n,a(n),i
call qsortint(iwork,n,a)
do i=1,n
iwork(n+i)=a(i)
enddo
do i=1,n
a(i)=iwork(n+iwork(i))
enddo
C
return
end
C
************************************************************************
SUBROUTINE QSORTINT(ORD,N,A)
************************************************************************
* SORTS THE ARRAY A(I),I=1,2,...,N BY PUTTING THE
* ASCENDING ORDER VECTOR IN ORD. THAT IS ASCENDING ORDERED A
* IS A(ORD(I)),I=1,2,...,N; DESCENDING ORDER A IS A(ORD(N-I+1)),
* I=1,2,...,N . THIS SORT RUNS IN TIME PROPORTIONAL TO N LOG N .
************************************************************************
IMPLICIT INTEGER (A-Z)
C
DIMENSION ORD(N),POPLST(2,20)
integer X,XX,Z,ZZ,Y
C
C TO SORT DIFFERENT INPUT TYPES, CHANGE THE FOLLOWING
C SPECIFICATION STATEMENTS; FOR EXAMPLE, FOR FORTRAN CHARACTER
C USE THE FOLLOWING: CHARACTER *(*) A(N)
C
integer A(N)
C
NDEEP=0
U1=N
L1=1
DO 1 I=1,N
1 ORD(I)=I
2 IF (U1.LE.L1) RETURN
C
3 L=L1
U=U1
C
C PART
C
4 P=L
Q=U
C FOR CHARACTER SORTS, THE FOLLOWING 3 STATEMENTS WOULD BECOME
C X = ORD(P)
C Z = ORD(Q)
C IF (A(X) .LE. A(Z)) GO TO 2
C
C WHERE "CLE" IS A LOGICAL FUNCTION WHICH RETURNS "TRUE" IF THE
C FIRST ARGUMENT IS LESS THAN OR EQUAL TO THE SECOND, BASED ON "LEN"
C CHARACTERS.
C
X=A(ORD(P))
Z=A(ORD(Q))
IF (X.LE.Z) GO TO 5
Y=X
X=Z
Z=Y
YP=ORD(P)
ORD(P)=ORD(Q)
ORD(Q)=YP
5 IF (U-L.LE.1) GO TO 15
XX=X
IX=P
ZZ=Z
IZ=Q
C
C LEFT
C
6 P=P+1
IF (P.GE.Q) GO TO 7
X=A(ORD(P))
IF (X.GE.XX) GO TO 8
GO TO 6
7 P=Q-1
GO TO 13
C
C RIGHT
C
8 Q=Q-1
IF (Q.LE.P) GO TO 9
Z=A(ORD(Q))
IF (Z.LE.ZZ) GO TO 10
GO TO 8
9 Q=P
P=P-1
Z=X
X=A(ORD(P))
C
C DIST
C
10 IF (X.LE.Z) GO TO 11
Y=X
X=Z
Z=Y
IP=ORD(P)
ORD(P)=ORD(Q)
ORD(Q)=IP
11 IF (X.LE.XX) GO TO 12
XX=X
IX=P
12 IF (Z.GE.ZZ) GO TO 6
ZZ=Z
IZ=Q
GO TO 6
C
C OUT
C
13 CONTINUE
IF (.NOT.(P.NE.IX.AND.X.NE.XX)) GO TO 14
IP=ORD(P)
ORD(P)=ORD(IX)
ORD(IX)=IP
14 CONTINUE
IF (.NOT.(Q.NE.IZ.AND.Z.NE.ZZ)) GO TO 15
IQ=ORD(Q)
ORD(Q)=ORD(IZ)
ORD(IZ)=IQ
15 CONTINUE
IF (U-Q.LE.P-L) GO TO 16
L1=L
U1=P-1
L=Q+1
GO TO 17
16 U1=U
L1=Q+1
U=P-1
17 CONTINUE
IF (U1.LE.L1) GO TO 18
C
C START RECURSIVE CALL
C
NDEEP=NDEEP+1
POPLST(1,NDEEP)=U
POPLST(2,NDEEP)=L
GO TO 3
18 IF (U.GT.L) GO TO 4
C
C POP BACK UP IN THE RECURSION LIST
C
IF (NDEEP.EQ.0) GO TO 2
U=POPLST(1,NDEEP)
L=POPLST(2,NDEEP)
NDEEP=NDEEP-1
GO TO 18
C
C END SORT
C END QSORT
C
END
C
************************************************************************
subroutine removed3(calc1)
************************************************************************
* Removes the -D3 postfix if it is specified in the calc or dft keywords
************************************************************************
implicit none
character*1 calc1(*)
integer i
c
do i=1,30
if(calc1(i ).eq.'-'.and.
$ calc1(i+1).eq.'d'.and.
$ calc1(i+2).eq.'3') then
calc1(i:i+2)=' '
exit
endif
enddo
c
return
end
C
************************************************************************
subroutine get_subsyscoord(natoms,atsymbol,coord,atnum,ind)
************************************************************************
* Write XYZ of the embedded subsystem
************************************************************************
#include "MRCCCOMMON"
integer, intent(in) :: natoms
integer, intent(in) :: atnum(*)
character(len=2), intent(in) :: atsymbol(natoms)
double precision, intent(in) :: coord(3,natoms)
integer, intent(inout) :: ind(natoms)
integer :: iatom,j,m
character(len=10) :: cscr
ind = 0
call embedat(natoms,minpfile,ind,iout,'embed ',5)
C Write COORD_SUBSYSA.xyz file
open(moldenfile,file='COORD_SUBSYSA.xyz')
m = 0
do iatom = 1 , natoms
if(atnum(iatom).gt.0.and.ind(iatom).eq.1) m=m+1
enddo
write(cscr,"(i10)") m
write(moldenfile,"(a10)") adjustl(cscr)
write(moldenfile,*)
do iatom = 1 , natoms
if(atnum(iatom).gt.0.and.ind(iatom).eq.1)
$write(moldenfile,"(a2,3f20.10)")
$adjustl(atsymbol(iatom)),(angtobohr*coord(j,iatom),j=1,3)
enddo
write(moldenfile,*)
close(moldenfile)
C
end subroutine get_subsyscoord
C
************************************************************************
character(len=16) function get_moselectalg( minpfile )
$ result( mo_select_alg )
************************************************************************
* Returns the MO selection scheme in the case of embedding calculations
************************************************************************
implicit none
integer, intent(in) :: minpfile
character(len=8) :: embed
character(len=32) :: dft
character(len=512) :: read_buffer
integer :: i,j
integer :: istat, channelopen
logical :: foundopen
istat = 0
read_buffer = ''
mo_select_alg = ' '
inquire(file='MINP',opened=foundopen,number=channelopen)
if( .not. foundopen .or. channelopen .ne. minpfile)
$ open(minpfile,file='MINP',position='REWIND',action='READ',
$ status='OLD',form='FORMATTED',iostat=istat)
call getkeym('embed',5,embed,8)
! atom list
read(minpfile,'(a)',iostat=istat) read_buffer
! low-level method
if(istat.eq.0) read(minpfile,'(a)',iostat=istat) read_buffer
dft = trim(adjustl( read_buffer ) )
if(trim(dft).eq.'user'.or.trim(dft).eq.'userd') then
read(minpfile,'(a)',iostat=istat) read_buffer
read( read_buffer , * ) j
do i=1,j
read(minpfile,'(a)',iostat=istat) read_buffer
enddo
if(trim(dft).eq.'userd') then
read(minpfile,'(a)',iostat=istat) read_buffer
read( read_buffer , * ) j
do i=1,j
read(minpfile,'(a)',iostat=istat) read_buffer
enddo
endif
if(istat.ne.0) then
write(6,'(a)') ' Error during the read of the
$user-defined DFT functional parameters'
call mrccend(1)
endif
endif
! Read the active MO specifier line
if(istat.eq.0) read(minpfile,'(a)',iostat=istat)
$ read_buffer
if( .not. foundopen .or. channelopen .ne. minpfile)
$close(minpfile)
! first character a number
i = index( read_buffer , ' ')
read_buffer( 1 : i ) = ' '
read_buffer = adjustl( read_buffer )
! second character: the method
i = index( read_buffer ,' ')
mo_select_alg = 'default'
if( i .ne. 1 ) mo_select_alg = read_buffer( 1 : i )
end function get_moselectalg
C
************************************************************************
subroutine get_orbloce_special(orblocc,orbloco,orblocv)
************************************************************************
#include "MRCCCOMMON"
integer :: istat, channelopen
character(len=16) :: cscr16
character(len=16) :: orblocc
character(len=16) :: orbloco
character(len=16) :: orblocv
logical :: foundopen
istat = 0
inquire(file='MINP',opened=foundopen,number=channelopen)
if(.not. foundopen .or. channelopen .ne. minpfile) then
open(unit=minpfile,file='MINP',status='OLD',iostat=istat,
$ action='READ',form='FORMATTED')
endif
rewind(minpfile)
call getkeym('orbloce',7,cscr16,16)
if(istat.eq.0) read(minpfile,'(a)',iostat=istat) orblocc
if(istat.eq.0) read(minpfile,'(a)',iostat=istat) orbloco
if(istat.eq.0) read(minpfile,'(a)',iostat=istat) orblocv
if(.not. foundopen .or. channelopen .ne. minpfile) then
close(minpfile)
endif
if( istat .ne. 0 ) then
write(iout,'(a)')
$'Cannot read orblocc/orbloco/orblocv for orbloce=special'
call mrccend(1)
endif
end subroutine get_orbloce_special
C
************************************************************************
subroutine separate_string(string,stringlen,separator,s1,s1len,
$s2,s2len)
************************************************************************
* Separates 'string' into string 's1' and 's2' connected by the separator
************************************************************************
implicit none
integer :: stringlen,i,j,s1len,s2len
character(len=1) :: line(stringlen),line2(stringlen)
character(len=1) :: s1(s1len),s2(s2len)
character(len=1) :: separator
character(len=1) :: string(stringlen)
logical llg
call lowercase(string,string,stringlen)
line(:)=string(:)
line2=' '
s1=' '
s2=' '
i=1
llg=.false.
do while(line(i).ne.' '.and.i.lt.stringlen)
if(line(i).eq.separator) then
line2(1:i-1)=line(1:i-1)
s1(1:s1len)=line2(1:s1len)
line2=' '
j=1
do while(line(i+j).ne.' '.and.j.ne.s2len.and.(i+j).ne.stringlen)
s2(j)=line(i+j)
j=j+1
enddo
llg=.true.
endif
i=i+1
enddo
if(.not.llg) then
s1(1:s1len)=string(1:s1len)
s2=' '
endif
return
end
C
************************************************************************
subroutine write_mop(qmmm,dens,isp,minpfile,natoms,natwdummy,
$ncent,coord,atsymbol,atchg,method,charge,mult,iout,
$usetemp,verbosity,uval,pcm,pcmkeystring,npcmdefkeys,
$oniomtype,ccoord,ncharge,firstcall,moldenval)
************************************************************************
* write MOPAC input file
************************************************************************
implicit none
integer ncent,natoms,isp,i,j,iout,minpfile,iatoms,k,ncharge
integer charge,mult,natwdummy,ierr,npcmdefkeys
character*1 line1(512)
character*2 atsymbol(natwdummy),cscr2
character*3 moldenval
character*4 verbosity,dens,uval,oniomtype
character*8 qmmm,cscr8
character*16 cscr16,inpfile
character*32 method,pcm
character*512 line,pcmkeystring
real*8 coord(3,ncent),atchg(ncent),ccoord(3,ncent),rscr
logical usetemp,llg,firstcall
equivalence(line,line1)
line=' '
write(iout,'(a)') ' Writing mop file ...'
write(cscr2,'(i2)') isp
! Search for template file
inpfile='MINP.'//trim(adjustl(cscr2))
cscr16='MINP.'//trim(adjustl(cscr2))//'.tpl'
if(usetemp) then
inquire(file=cscr16,exist=llg)
if(llg) then
write(iout,
$"( ' Template file is found. Using file: ',a12)") cscr16
open(minpfile, file=cscr16, iostat=ierr)
rewind(minpfile)
read(minpfile,'(a)') line
close(minpfile)
else
write(iout,"( ' Template file is not found.')")
endif
endif
line=trim(adjustl(line))
! set MS
rscr=(mult-1d0)/2.0d0
write(cscr16,'(f6.1)') rscr
line=trim(adjustl(line))//' MS='//trim(adjustl(cscr16))
if(mult.ne.1) line=trim(adjustl(line))//' UHF '
! set Charge
write(cscr16,'(i5)') charge
line=trim(adjustl(line))//' CHARGE='//trim(adjustl(cscr16))
! set Bohr/Angstrom
if(uval.eq.'bohr') line=trim(adjustl(line))//' A0 '
! set verbosity
read(verbosity,'(i4)') i
if(i.gt.2) line=trim(adjustl(line))//' '//'PL'
! set bond order
if(isp.eq.1) line=trim(adjustl(line))//' BONDS'
! Set gradient and density-reading related keywords
read(dens,'(i4)') i
! Set gradient and density matrix writing keywords
if(i.gt.1) line=trim(adjustl(line))//' GRADIENTS DENOUT'
! Set density matrix reading keyword (use it as guess)
if(.not.firstcall) line=trim(adjustl(line))//' OLDENS '
! Print MOs in the unformatted, Molden-readable format
! (the formatted command, GRAPHF, is not working... )
if(moldenval.ne.'off') line=trim(adjustl(line))//' GRAPH '
!!! Most important keywords for the interface !!!
! Request a single SCF
line=trim(adjustl(line))//' 1SCF '
! Request an auxilirary file with increased printing precision
line=trim(adjustl(line))//' AUX(PRECISION=12,XP,XS,XW) '
! Disable the reorientation of the system
line=trim(adjustl(line))//' NOREOR '
! Increase the SCF criteria (do we need this other than developement?)
line=trim(adjustl(line))//' SCFCRT=1.D-8 '
! Add the calculation type
line=trim(adjustl(method))//' '//trim(adjustl(line))
! Add solvent-specification
if(npcmdefkeys.gt.0) then
line=trim(adjustl(line))//' '//trim(adjustl(pcmkeystring))
endif
open(minpfile, file=inpfile, iostat=ierr, status='replace')
write(minpfile,'(a)') trim(adjustl(line))
write(minpfile,*)
if(usetemp) then
write(minpfile,'(a)')
$'! User defined keywords from: template file+MINP'
else
write(minpfile,'(a)')'! User defined keywords from: MINP'
endif
do iatoms=1,natoms
write(minpfile,'(a2,3f25.16)')
$atsymbol(iatoms),coord(1:3,iatoms)
enddo
close(minpfile)
return
end
C
************************************************************************
subroutine execute_mopac(cisp,iprog,iout,istat,cput,walt)
************************************************************************
* execute MOPAC
************************************************************************
implicit none
#if !defined (gfortran) && !defined (G95)
integer*4 system
external system
#endif
integer iout,istat
character*2 cisp
character*16 iprog
character*64 cscr64
real*8 cput,walt
write(iout,'(2a)')
$' MOP file is written into: MINP.'//trim(adjustl(cisp))
cscr64=trim(adjustl(iprog))//' MINP.'//trim(adjustl(cisp))
write(iout,'(2a)')
$' Executing command line: '//trim(adjustl(cscr64))
cput=0.0d0
walt=0.0d0
call mtime(cput,walt)
istat=system(trim(adjustl(cscr64)))
cput=1.0d0
walt=1.0d0
call mtime(cput,walt)
return
end
************************************************************************
subroutine read_mopac_aux(cfilenum,scrfile1,nval,rscrmat,rtype,
$skey,skeylen)
************************************************************************
* read MOPAC's aux file
************************************************************************
implicit none
integer :: ios,scrfile1,skeylen,i,n,nval,nrow
real*8 :: rscrmat(nval),rscrr
character(len=2) :: cfilenum
character(len=6) :: rtype
character(len=1) :: bline(64),linev,cscr1
character(len=64) :: read_buffer,cscr64
character(len=64) :: skey
equivalence(read_buffer,bline)
equivalence(cscr64,linev)
open(scrfile1,file='MINP.'//trim(adjustl(cfilenum))//'.aux',
$ form='formatted')
rewind(scrfile1)
cscr64=skey
linev(skeylen:64)=' '
do
read(scrfile1,'(a)',iostat=ios) read_buffer
! End of file; data not found
bline(skeylen:64)=' '
if(ios.lt.0) exit
if (read_buffer.eq.cscr64) then
if(rtype.eq.'block ') then
read(scrfile1,*) cscr1
if(cscr1.ne.'#') backspace(scrfile1)
! Count rows
rscrr=nval/10d0
nrow=nint(rscrr)
if(nval.gt.0.and.nrow.eq.0) nrow=1
! Read from aux file
do i=1,nrow
if(i.ne.nrow) then
read(scrfile1,*) rscrmat((i-1)*10+1:10*i)
else
read(scrfile1,*) rscrmat((i-1)*10+1:nval)
endif
enddo
else if(rtype.eq.'horizo') then
backspace(scrfile1)
read(scrfile1,'(a)') read_buffer
i=1
do while(bline(i).ne.'='.and.i.le.64)
bline(i)=' '
i=i+1
enddo
bline(i)=' '
read(read_buffer,*) rscrmat(1:nval)
endif
endif
enddo
close(scrfile1)
return
end
************************************************************************
subroutine check_sqm_prog(sqmprog,exe,iout,userset)
************************************************************************
* Check the availability of sqm programs
************************************************************************
implicit none
character*16 sqmprog,exe,extprog(2),extprogexe(2)
character*32 string
integer iout,istat,nprog,iprog,foundprog(2)
logical lfound,userset
#if !defined (gfortran) && !defined (G95)
integer*4 system
external system
#endif
lfound=.false.
nprog=2
foundprog=0
extprog(1)='xtb '
extprog(2)='mopac2016 '
extprogexe(1)='xtb '
extprogexe(2)='MOPAC2016.exe '
write(iout,*)
write(iout,'(a)')
$' Checking if a semi-empirical program is '//
$'available...'
do iprog=1,nprog
if(sqmprog.ne.' '.and.
$ sqmprog.ne.extprog(iprog)) then
cycle
endif
write(iout,"(' Searching for sqm program: ',a)")
$ trim(adjustl(extprog(iprog)))
string='which '//trim(adjustl(extprogexe(iprog)))
istat=system(trim(adjustl(string)))
foundprog(iprog)=1
if(istat.eq.0) then
lfound=.true.
write(iout,"(' 'a,' is found!')")
$ trim(adjustl(extprogexe(iprog)))
else
lfound=.false.
endif
enddo
if(.not.lfound) then
write(iout,"(' No executables are found.')")
call mrccend(1)
else
do iprog=1,nprog
if((sqmprog.eq.' '.and.foundprog(iprog).eq.1)
$.or.
$ (sqmprog.eq.extprog(iprog) .and.foundprog(iprog).eq.1))
$then
exe=extprogexe(iprog)
exit
else if(sqmprog.eq.extprog(iprog) .and.foundprog(iprog).eq.0)
$ then
write(iout,"(' The executable for the requested program ('a,
$') is not found.')") trim(adjustl(sqmprog))
call mrccend(1)
endif
enddo
write(iout,*)
write(iout,"(' ',a,' will be used as external sqm program.')")
$ trim(adjustl(exe))
if(.not.userset) then
write(iout,'(a)') ' (by default) '
else
write(iout,'(a)') ' (by user request) '
endif
write(iout,*)
endif
return
end
************************************************************************
subroutine write_xtb_files(qmmm,dens,isp,minpfile,natoms,
$natwdummy,ncent,coord,atsymbol,atchg,method,charge,mult,
$iout,usetemp,verbosity,uval,pcm,pcmkeystring,npcmdefkeys,
$oniomtype,ccoord,ncharge,firstcall,moldenval,xtbstring,
$printcharges,initguess,etemp)
************************************************************************
* write input files for the xtb program
************************************************************************
implicit none
integer ncent,natoms,isp,i,j,iout,minpfile,iatoms,k,ncharge
integer charge,mult,natwdummy,ierr,npcmdefkeys
character*1 line1(512)
character*2 atsymbol(natwdummy),cscr2
character*3 moldenval
character*4 verbosity,dens,uval,oniomtype
character*8 qmmm,cscr8,etemp
character*16 cscr16,inpfile,pcfile,coordfile,coordfile2
character*16 tempfile,outfile
character*32 method,pcm,cscr32
character*512 line,pcmkeystring
character*128 xtbstring
real*8 coord(3,ncent),atchg(ncent),ccoord(3,ncent),rscr
logical usetemp,llg,firstcall,printcharges,initguess
equivalence(line,line1)
line=' '
llg=.false.
write(iout,'(a)') ' Writing files for xtb program ...'
write(cscr2,'(i2)') isp
! Search for template file
inpfile='MINP.'//trim(adjustl(cscr2))//'.xtbinp'
coordfile='MINP.'//trim(adjustl(cscr2))//'.xyz'
coordfile2='MINP.'//trim(adjustl(cscr2))//'.coord'
pcfile='MINP.'//trim(adjustl(cscr2))//'.pc'
tempfile='MINP.'//trim(adjustl(cscr2))//'.tpl'
! Read the first line of the template file to extend the command line
if(usetemp) then
inquire(file=tempfile,exist=llg)
if(llg) then
write(iout,
$"( ' Template file is found. Using file: ',a12)") tempfile
open(minpfile, file=tempfile, iostat=ierr)
rewind(minpfile)
read(minpfile,'(a)') line
close(minpfile)
else
write(iout,"( ' Template file is not found.')")
endif
endif
xtbstring='--input '//trim(adjustl(inpfile))//' --grad '//
$'--namespace MINP.'//trim(adjustl(cscr2))//' '//
$trim(adjustl(line))
if(initguess) then
xtbstring=trim(adjustl(xtbstring))//' --acc 0.0001'
endif
xtbstring=trim(adjustl(xtbstring))//' --etemp '//
$trim(adjustl(etemp))
open(minpfile,file=inpfile,form='formatted')
! set spin
c rscr=(mult-1d0)/2.0d0
c write(cscr16,'(f6.1)') rscr
i=(mult-1) ! xtb requires the number of unpaired electrons
write(cscr16,'(i4)') i
write(minpfile,'(2a)') '$spin '//trim(adjustl(cscr16))
! set Charge
write(cscr16,'(i5)') charge
write(minpfile,'(2a)') '$chrg '//trim(adjustl(cscr16))
! set geometry units - geometry input dependent??
if(uval.eq.'bohr') then
write(minpfile,'(a)') '$coord bohr'
else
write(minpfile,'(a)') '$coord angs'
endif
write(minpfile,'(a)') '$write'
! set bond order
if(isp.eq.1) write(minpfile,'(a)') ' wiberg=1'
! Do not print (mulliken?) charges
if(printcharges) then
write(minpfile,'(a)') ' charges=1'
else
write(minpfile,'(a)') ' charges=0'
endif
! set molden output
if(moldenval.ne.'off') write(minpfile,'(a)') ' mos=1'
write(minpfile,'(a)') '$'
! Set gradient and density-reading related keywords
read(dens,'(i4)') i
! Set gradient and density matrix writing keywords
if(i.gt.1) xtbstring=trim(adjustl(xtbstring))//' --grad'
! Add the calculation type
if(method.eq.'gfn0-xtb'.or.
$ method.eq.'gfn1-xtb'.or.
$ method.eq.'gfn2-xtb') then
write(minpfile,'(a)') '$gfn'
if(method.eq.'gfn0-xtb') write(minpfile,'(a)') 'method=0'
if(method.eq.'gfn1-xtb') write(minpfile,'(a)') 'method=1'
if(method.eq.'gfn2-xtb') write(minpfile,'(a)') 'method=2'
write(minpfile,'(a)') '$'
else if(method.eq.'gfn-ff') then
xtbstring=trim(adjustl(xtbstring))//' --gfnff'
endif
if(ncharge.gt.0) then
write(minpfile,'(a)') '$embedding'
write(minpfile,'(3a)')
$ 'input=MINP.'//trim(adjustl(cscr2))//'.pc'
write(minpfile,'(a)') '$'
endif
close(minpfile)
! Catanate the template file and the input file
if(usetemp.and.llg) then
call ishell('cat '//trim(adjustl(tempfile))//' >>'//
$trim(adjustl(inpfile)))
endif
open(minpfile, file=coordfile, iostat=ierr, status='replace')
write(cscr16,'(i16)') natoms
write(minpfile,'(a)') trim(adjustl(cscr16))
write(minpfile,*)
do iatoms=1,natoms
write(minpfile,'(a2,3f25.16)')
$atsymbol(iatoms),coord(1:3,iatoms)
enddo
write(minpfile,*)
close(minpfile)
!
open(minpfile, file=coordfile2, iostat=ierr, status='replace')
if(uval.eq.'angs') write(minpfile,'(a)') '$coord angs'
if(uval.eq.'bohr') write(minpfile,'(a)') '$coord bohr'
do iatoms=1,natoms
write(minpfile,'(3f25.16,2x,a2)')
$coord(1:3,iatoms),atsymbol(iatoms)
enddo
write(minpfile,'(a)') '$end'
close(minpfile)
! Write point charge file in the case of embedding
if(ncharge.gt.0) then
open(minpfile, file=pcfile, iostat=ierr, status='replace')
write(minpfile,'(100000000i12)') ncharge
do iatoms=1,ncharge
write(minpfile,'(4f25.16,2x,a)')
$atchg(iatoms),ccoord(1:3,iatoms),' 999'
enddo
close(minpfile)
endif
! Add pcm info
if(trim(adjustl(pcm)).ne.'off') then
if(npcmdefkeys.ne.0) then
do i=1,npcmdefkeys
if(i.eq.1) then
read(pcmkeystring,*) line
else
read(pcmkeystring,*) (line1(k),k=1,i-1),line
endif
call separate_string(line,512,'=',cscr32,32,cscr16,16)
if (trim(adjustl(cscr16)).eq.'alpb') then
xtbstring=trim(adjustl(xtbstring))//' --alpb '//pcm
else if(trim(adjustl(cscr16)).eq.'gbsa') then
xtbstring=trim(adjustl(xtbstring))//' --gbsa '//pcm
endif
enddo
else
xtbstring=trim(adjustl(xtbstring))//' --alpb '//pcm
endif
endif
! Add the coordinate file to the executing string
xtbstring='xtb '
$//trim(adjustl(coordfile2))//' '//trim(adjustl(xtbstring))
return
end
C
************************************************************************
subroutine execute_xtb(cisp,iout,xtbstring,ioutprint,istat,
$cput,walt)
************************************************************************
* execute MOPAC
************************************************************************
implicit none
integer iout,istat
#if !defined (gfortran) && !defined (G95)
integer*4 system
external system
#endif
character*2 cisp
character*16 outfile
character*128 xtbstring
real*8 cput,walt
logical ioutprint
outfile='MINP.'//trim(adjustl(cisp))//'.xtbout'
write(iout,'(2a)')
$' xTB xcontrol file is written into: MINP.'//trim(adjustl(cisp))
$//'.xtbinp'
write(iout,'(3a)')
$' xTB coord file is written into : '//
$'MINP.'//trim(adjustl(cisp))//'.xyz'
if(.not.ioutprint) then
xtbstring=trim(adjustl(xtbstring))//' > '
$//trim(adjustl(outfile))
write(iout,'(3a)')
$' xTB output file is written into : '//
$'MINP.'//trim(adjustl(cisp))//'.xtbout'
endif
write(iout,'(a)') ' Executing command line: '
write(iout,'(x,a)') trim(adjustl(xtbstring))
cput=0.0d0
walt=0.0d0
call mtime(cput,walt)
istat=system(trim(adjustl(xtbstring)))
cput=1.0d0
walt=1.0d0
call mtime(cput,walt)
return
end
C
logical function mopac_was_used() result( used )
#include "MRCCCOMMON"
character(len=1024) :: read_buffer
character(len=8) :: prop
character(len=16) :: method
integer :: istat
used = .false.
open(unit=ifcfile,file='iface',status='OLD',
$ action='READ',position='REWIND',iostat=istat)
read(ifcfile,'(a)',iostat=istat) read_buffer
do while( istat .eq. 0 .and. .not. used )
read(ifcfile,'(a)',iostat=istat) read_buffer
read(read_buffer,*) prop,method
used = trim(adjustl(method)) .eq. 'MOPAC-SCF'
enddo
close(ifcfile)
end function mopac_was_used
************************************************************************
subroutine getatnum(atsymbol,atnum)
************************************************************************
* Get atomic number and atomic symbol
************************************************************************
#include "MRCCCOMMON"
integer atnum
character*2 atsymbol,atsym(0:natmax),ch1,ch2
data atsym /'Bf',
$ 'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne',
$ 'Na','Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca',
$ 'Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn',
$ 'Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y ','Zr',
$ 'Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn',
$ 'Sb','Te','I ','Xe','Cs','Ba','La','Ce','Pr','Nd',
$ 'Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb',
$ 'Lu','Hf','Ta','W ','Re','Os','Ir','Pt','Au','Hg',
$ 'Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th',
$ 'Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es','Fm',
$ 'Md','No','Lr','Rf','Db','Sg','Bh','Hs','Mt','Ds',
$ 'Rg','Cn','Nh','Fl','Mc','Lv','Ts','Og'/
C
if(atnum.eq.0) then
call lowercase(atsymbol,ch1,2)
atnum=-1
ch2=' '
do while(ch1.ne.ch2)
atnum=atnum+1
if(atnum.gt.natmax) then
call uppercase(atsymbol,atsymbol,1)
write(iout,*) 'Unknown atomic symbol '//trim(atsymbol)//'!'
call mrccend(1)
endif
call lowercase(atsym(atnum),ch2,2)
enddo
endif
atsymbol=atsym(atnum)
C
return
end
C
***********************************************************************
logical function dof12(calc)
***********************************************************************
* Shall we do F12 calculation?
***********************************************************************
implicit none
character*1 calc(16)
integer i
C
dof12=.false.
do i=1,13
dof12=dof12.or.((calc(i ).eq.'-'.or.calc(i ).eq.'(').and.
$ calc(i+1).eq.'f'.and.
$ calc(i+2).eq.'1'.and.
$ calc(i+3).eq.'2')
enddo
C
return
end
C
************************************************************************
subroutine orth(nb,mb,s,c,w)
************************************************************************
* Symmetric orthogonalization
************************************************************************
implicit none
integer nb,mb,i,j
real*8 s(nb,nb),c(nb,mb),w(*)
C
call dsymm('l','u',nb,mb,1.d0,s,nb,c,nb,0.d0,w(nb*nb+1),nb)
call dgemm('t','n',mb,mb,nb,1.d0,c,nb,w(nb*nb+1),nb,0.d0,w,mb)
call invsqrt(w,mb,w(mb*mb+1),6,i,0.d0,1)
if(i.ne.0) call mrccend(1)
call dsymm('r','u',nb,mb,1.d0,w,mb,c,nb,0.d0,w(mb*mb+1),nb)
call dcopy(nb*mb,w(mb*mb+1),1,c,1)
C
return
end
C
************************************************************************
subroutine smearoccup(nmos,eigenvalue,nelectr,occup,etemp)
************************************************************************
* Calculate the occupation of the orbitals if the Fermi-smearing is utilized
************************************************************************
#include "MRCCCOMMON"
integer nmos,nelectr,p,nsingle
real*8 fermilevel,betafact,etemp,sumocc,sumscal,sumfod
real*8 occup(nmos),eigenvalue(nmos)
sumfod=0.0d0
nsingle=0
sumocc=0.0d0
occup(1:nmos)=0.0d0
! Calculate fractional occupations
fermilevel=( eigenvalue(nelectr)+eigenvalue(nelectr+1) )/2.d0
betafact=1.0d0/((boltz/jtoeh*1d-5)*etemp)
do p=1,nmos
occup(p)=1.0d0/
$(1.0d0+exp( (eigenvalue(p)-fermilevel)*betafact ) )
sumocc=sumocc+occup(p)
if(abs(1.0d0-occup(p)).lt.1d-9) then
nsingle=nsingle+1
else
sumfod=sumfod+occup(p)
endif
enddo
! Rescale the occupation number of the partially occupied orbitals
sumscal=nsingle
if(sumfod.gt.1d-9) then
do p=nsingle+1,nmos
occup(p)=dble(nelectr-nsingle)/sumfod*occup(p)
sumscal=sumscal+occup(p)
enddo
endif
c write(iout,'(a,E15.10)') ' Particle number residue: ',
c $abs( sumscal-dble(nelectr) )
return
end
C
************************************************************************
subroutine gsch(o, x, y, ldo, n, m)
************************************************************************
* gram-schmidt ortogonalization
* input:
* o(ldo, m) a set of orthonormal vectors
* y: temporary vector of length ldo
* x: new vector
* output:
* x: ortogonalized vector
************************************************************************
implicit none
integer ldo, n, m, i
double precision o(ldo, *), x(*), norm, y(*), tmp, oldnorm, eta
parameter(eta = 0.707106781d0)
double precision dnrm2, ddot
if(m .ge. n) then
call dlaset('f', 1, n, 0.0d0, 0.0d0, x, 1)
return
endif
oldnorm = dnrm2(n, x, 1)
call dgemv('t', n, m, 1.0d0, o, ldo, x, 1, 0.0d0, y, 1)
call dgemv('n', n, m, -1.0d0, o, ldo, y, 1, 1.0d0, x, 1)
norm = dnrm2(n, x, 1)
do while(oldnorm*eta .gt. norm .or.
& (oldnorm*eta .gt. norm .and. norm .gt. 1.0d-16))
call dgemv('t', n, m, 1.0d0, o, ldo, x, 1, 0.0d0, y, 1)
call dgemv('n', n, m, -1.0d0, o, ldo, y, 1, 1.0d0, x, 1)
oldnorm = norm
norm = dnrm2(n, x, 1)
enddo
norm = dnrm2(n, x, 1)
if(norm .gt. 1.0d-16) then
call dscal(n, 1.0d0/norm, x, 1)
else
call dlaset('f', 1, n, 0.0d0, 0.0d0, x, 1)
endif
end subroutine
************************************************************************
subroutine write_dat_file_for_amber
$(finalener,selfener,idens,natoms,ncent,grads,efield,
$dipx,dipy,dipz)
************************************************************************
* Write a formatted data file for the Amber interface
************************************************************************
#include "MRCCCOMMON"
integer :: datoms,natoms,ncent,idens
real*8 :: grads(3,*)
real*8 :: efield(3,*)
real*8 :: selfener
real*8 :: finalener
real*8 :: dipx,dipy,dipz
C
open(scrfile1,file='mrcc_job.dat',form='formatted')
if(idens.gt.1) then
write(scrfile1,'(a)')
write(scrfile1,'(a)') 'Forces on QM atoms [Hartree/Bohr]'
do datoms=1,natoms
write(scrfile1,"(3e23.15)")
$ grads(1,datoms),grads(2,datoms),grads(3,datoms)
enddo
C Electric field
write(scrfile1,'(a)')
write(scrfile1,'(a)') 'Electric field at MM atoms'
do datoms=natoms+1,ncent
write(scrfile1,"(3e23.15)")
$ efield(1,datoms),efield(2,datoms),efield(3,datoms)
enddo
C Dipole moment
write(scrfile1,'(a)')
write(scrfile1,'(a)')
$'Magnitude of dipole moment vector [Debye]'
write(scrfile1,'(e23.15)')
$echesu*angtobohr*dsqrt(dipx**2+dipy**2+dipz**2)
write(scrfile1,'(a)')
$'Direction of dipole moment vector [Debye]'
write(scrfile1,'(3e23.15)')
$echesu*angtobohr*dipx,echesu*angtobohr*dipy,echesu*angtobohr*dipz
write(scrfile1,'(a)')
endif
write(scrfile1,'(a)')'Final energy of the QM region [AU]'
write(scrfile1,'(e23.15)') finalener
write(scrfile1,'(a)')'Self energy of the point charges [AU]'
write(scrfile1,'(e23.15)') selfener
close(scrfile1)
end subroutine write_dat_file_for_amber
C
************************************************************************
subroutine readh(r8heap,ifile,symtra,nbasis,offset,lsa)
************************************************************************
* Read core Hamiltonian
************************************************************************
implicit none
integer ifile,nbasis,offset
real*8 r8heap(*),symtra
logical lsa
C
open(ifile,file='OEINT',form='unformatted')
read(ifile)
call readone(r8heap,ifile,symtra,nbasis,r8heap(nbasis**2+1),
$offset,lsa)
close(ifile)
C
return
end
C
************************************************************************
subroutine readone(mat,ifile,symtra,nbasis,r8heap,offset,lsa)
************************************************************************
* Read and compress one-electron integrals
************************************************************************
implicit none
integer ifile,nbasis,offset,ii
real*8 mat(nbasis,nbasis),symtra(nbasis,nbasis),r8heap(*)
logical lsa
C
if(lsa) then
ii=nbasis*nbasis+1
call roeint(r8heap(ii),r8heap(ii),r8heap,ifile,nbasis)
call mxts(r8heap,symtra,r8heap(ii),nbasis,mat,offset,lsa)
else
call roeint(r8heap,r8heap,mat,ifile,nbasis)
endif
C
return
end
C
************************************************************************
subroutine mxts(mat1,symtra,r8heap,nbasis,mat2,offset,lsa)
************************************************************************
* Transform one-electron quantities to symmetry basis
************************************************************************
implicit none
integer nbasis,offset
real*8 mat1(nbasis,nbasis),mat2(nbasis,nbasis)
real*8 symtra(nbasis,nbasis),r8heap(*)
logical lsa
C
if(.not.lsa) return
call mx_basis_tr(mat1,symtra,r8heap,'TS',nbasis)
call mx_symm_compr(mat1,mat2,offset)
C
return
end
C
************************************************************************
subroutine mxto(mat1,symtra,r8heap,nbasis,mat2,offset,lsa)
************************************************************************
* Transform one-electron quantities to symmetry basis
************************************************************************
implicit none
integer nbasis,offset
real*8 mat1(nbasis,nbasis),mat2(nbasis,nbasis)
real*8 symtra(nbasis,nbasis),r8heap(*)
logical lsa
C
if(.not.lsa) return
call mx_symm_extr(mat1,mat2,offset)
call mx_basis_tr(mat1,symtra,r8heap,'TO',nbasis)
C
return
end
C
************************************************************************
subroutine mx_basis_tr(mx,symtra,r8heap,mode,lsize)
************************************************************************
* subroutine to symmetrise matrices with transformation
* matrix from symtra file (running in destructive mode)
************************************************************************
implicit none
c arguments
real*8 mx(lsize**2) !non compressed input/output matrix
real*8 symtra(lsize**2) !symmetry transformation matrix
real*8 r8heap(*) !workspace
character*2 mode
integer lsize !linear size of matrices
c variables
character*1 lmx,rmx
c program
select case (mode)
case ('ts','TS','Ts','tS')
lmx='t'
rmx='n'
case ('to','TO','To','tO')
lmx='n'
rmx='t'
case default
write(*,*) ' @mx_basis_tr: ERROR'
write(*,*) ' unknown transformation mode'
write(*,*) ' last argument of subroutine'
return
end select
c multiply from right
call dgemm(
& 'n',
& rmx,
& lsize,
& lsize,
& lsize,
& 1.d0,
& mx,
& lsize,
& symtra,
& lsize,
& 0.d0,
& r8heap,
& lsize)
c multiply from left
call dgemm(
& lmx,
& 'n',
& lsize,
& lsize,
& lsize,
& 1.d0,
& symtra,
& lsize,
& r8heap,
& lsize,
& 0.d0,
& mx,
& lsize)
c
return
end
************************************************************************
subroutine mx_symm_extr(full_mx,compr_mx,offset)
************************************************************************
* subroutine to extract a matrix from blockdiagonal storage format
************************************************************************
#include "MRCCCOMMON"
#include "SCFCOMMON"
c arguments
real*8 full_mx(nbasis,nbasis)
real*8 compr_mx(sqrsize)
integer offset(*)
c scalar variables
integer pivot
integer ii,jj,kk
c program
call dfillzero(full_mx,nbasis**2)
pivot=1
do kk=1,nir
if(nfunc(kk).eq.0) cycle
do jj=1,nfunc(kk)
do ii=1,nfunc(kk)
full_mx(offset(kk)-1+ii,offset(kk)-1+jj)=compr_mx(pivot)
pivot=pivot+1
enddo
enddo
enddo
c
return
end
************************************************************************
subroutine mx_symm_compr(full_mx,compr_mx,offset)
************************************************************************
* subroutine to comprass a matrix into blockdiagonal storage format
************************************************************************
#include "MRCCCOMMON"
#include "SCFCOMMON"
c arguments
real*8 full_mx(nbasis,nbasis)
real*8 compr_mx(sqrsize)
integer offset(*)
c scalar variables
integer pivot
integer ii,jj,kk
c program
pivot=1
do kk=1,nir
if(nfunc(kk).eq.0) cycle
do jj=1,nfunc(kk)
do ii=1,nfunc(kk)
compr_mx(pivot)=full_mx(offset(kk)-1+ii,offset(kk)-1+jj)
pivot=pivot+1
enddo
enddo
enddo
c
return
end
************************************************************************
subroutine nactorb(corb,n)
************************************************************************
* Calculates the number of active/doubly occupied MCSCF orbitals
************************************************************************
implicit none
integer n,m,i,ii
character(len=1) line2(16),line3(512)
character(len=16) line1
character(len=512) line
character(len=4) corb
equivalence(line,line3)
equivalence(line1,line2)
C
line3=' '
call getkey(corb,4,line3,512)
n=0
ii=0
line2=' '
do i=1,len(line)
select case(line3(i))
case ('0','1','2','3','4','5','6','7','8','9')
ii=ii+1
line2(ii)=line3(i)
case (',')
read(line1,*) m
n=n+m
ii=0
line2=' '
case default
exit
end select
enddo
read(line1,*) m
n=n+m
C
return
end
C
************************************************************************
module error_handler
************************************************************************
* Error handling for various run-time problems.
* It is intended to give more informative feedbacks to the user
* while keeping programming duties minimal.
* 2023-05-12 Bence Hégely
************************************************************************
save
integer :: out_channel
interface out_of_range_error
module procedure out_of_range_error_limited_int
module procedure out_of_range_error_unlimited_int
module procedure out_of_range_error_limited_real
module procedure out_of_range_error_unlimited_real
module procedure out_of_range_error_limited_mixed_lowreal
module procedure out_of_range_error_limited_mixed_lowint
end interface
interface illegal_value_error
module procedure illegal_value_error_int
module procedure illegal_value_error_char
end interface
!
contains
!
subroutine init_error_handler( iout )
out_channel = iout
end subroutine
!
subroutine failed_memop(varname,memoptype)
implicit none
integer :: memoptype
character(len=*) :: varname
if(memoptype.eq.1) then
write(out_channel,"(a,a,a)")
$" Memory allocation is failed for variable: '",
$trim(adjustl(varname)),"'"
else if(memoptype.eq.2) then
write(out_channel,"(a,a,a)")
$" Memory deallocation is failed for variable: '",
$trim(adjustl(varname)),"'"
endif
call mrccend(1)
end subroutine failed_memop
!
subroutine io_error(inpstring,source)
implicit none
character(len=*) :: inpstring
character(len=*) :: source
write(out_channel,'(a)') ' An error occured during IO operation!'
write(out_channel,'(5a)') ' (',inpstring,' @ ',source,')'
call mrccend(1)
end subroutine io_error
!
subroutine out_of_range_error_limited_int
$(message,upper,ubound_in,lower,lbound_in)
implicit none
character(len=*),intent(in) :: message
integer, intent(in) :: upper
integer, intent(in) :: lower
logical, intent(in) :: ubound_in
logical, intent(in) :: lbound_in
character(len=32) :: cupper
character(len=32) :: clower
character(len=1 ) :: c1u
character(len=1 ) :: c1l
write(out_channel,'(a)') ' Out of bounds error.'
write(out_channel,'(a)') message
write(out_channel,'(a)',advance='no') ' Legal range: '
if( lbound_in) c1l = '['
if(.not.lbound_in) c1l = '('
if( ubound_in) c1u = ']'
if(.not.ubound_in) c1u = ')'
write(cupper,*) upper
write(clower,*) lower
write(out_channel,*)
$ c1l,trim(adjustl(clower)),',',trim(adjustl(cupper)),c1u
call mrccend(1)
end subroutine out_of_range_error_limited_int
!
subroutine out_of_range_error_unlimited_int
$(message,bound_type,bound,bound_in)
implicit none
character(len=*),intent(in) :: message
character(len=1 ),intent(in) :: bound_type
integer, intent(in) :: bound
logical, intent(in) :: bound_in
character(len=32) :: cbound
character(len=1 ) :: c1b
write(out_channel,'(a)') ' Out of bounds error.'
write(out_channel,'(a)') message
write(out_channel,'(a)',advance='no') ' Legal range: '
if( bound_type .eq. 'u' ) then
if( bound_in) c1b = ']'
if(.not.bound_in) c1b = ')'
write(cbound,*) bound
write(out_channel,*) '(-INFINITY,',trim(adjustl(cbound)),c1b
else if( bound_type .eq. 'l' ) then
if( bound_in) c1b = '['
if(.not.bound_in) c1b = '('
write(cbound,*) bound
write(out_channel,*) c1b,trim(adjustl(cbound)),',+INFINITY)'
else
write(6,*) ' Cannot show bounds.'
endif
call mrccend(1)
end subroutine out_of_range_error_unlimited_int
!
subroutine out_of_range_error_limited_real
$(message,upper,ubound_in,lower,lbound_in)
implicit none
character(len=*), intent(in) :: message
double precision, intent(in) :: upper
double precision, intent(in) :: lower
logical, intent(in) :: ubound_in
logical, intent(in) :: lbound_in
character(len=32) :: cupper
character(len=32) :: clower
character(len=1 ) :: c1u
character(len=1 ) :: c1l
write(out_channel,'(a)') ' Out of bounds error.'
write(out_channel,'(a)') message
write(out_channel,'(a)',advance='no') ' Legal range: '
if( lbound_in) c1l = '['
if(.not.lbound_in) c1l = '('
if( ubound_in) c1u = ']'
if(.not.ubound_in) c1u = ')'
write(cupper,*) upper
write(clower,*) lower
write(out_channel,*)
$ c1l,trim(adjustl(clower)),',',trim(adjustl(cupper)),c1u
call mrccend(1)
end subroutine out_of_range_error_limited_real
!
subroutine out_of_range_error_unlimited_real
$(message,bound_type,bound,bound_in)
implicit none
character(len=*), intent(in) :: message
character(len=1 ), intent(in) :: bound_type
double precision , intent(in) :: bound
logical, intent(in) :: bound_in
character(len=32) :: cbound
character(len=1 ) :: c1b
write(out_channel,'(a)') ' Out of bounds error.'
write(out_channel,'(a)') message
write(out_channel,'(a)',advance='no') ' Legal range: '
if( bound_type .eq. 'u' ) then
if( bound_in) c1b = ']'
if(.not.bound_in) c1b = ')'
write(cbound,*) bound
write(out_channel,*) '(-INFINITY,',trim(adjustl(cbound)),c1b
else if( bound_type .eq. 'l' ) then
if( bound_in) c1b = '['
if(.not.bound_in) c1b = '('
write(cbound,*) bound
write(out_channel,*) c1b,trim(adjustl(cbound)),',+INFINITY)'
else
write(out_channel,*) ' Cannot show bounds.'
endif
call mrccend(1)
end subroutine out_of_range_error_unlimited_real
!
subroutine out_of_range_error_limited_mixed_lowreal
$(message,upper,ubound_in,lower,lbound_in)
implicit none
character(len=*), intent(in) :: message
integer , intent(in) :: upper
double precision, intent(in) :: lower
logical, intent(in) :: ubound_in
logical, intent(in) :: lbound_in
character(len=32) :: cupper
character(len=32) :: clower
character(len=1 ) :: c1u
character(len=1 ) :: c1l
write(out_channel,'(a)') ' Out of bounds error.'
write(out_channel,'(a)') message
write(out_channel,'(a)',advance='no') ' Legal range: '
if( lbound_in) c1l = '['
if(.not.lbound_in) c1l = '('
if( ubound_in) c1u = ']'
if(.not.ubound_in) c1u = ')'
write(cupper,*) upper
write(clower,*) lower
write(out_channel,*)
$ c1l,trim(adjustl(clower)),',',trim(adjustl(cupper)),c1u
call mrccend(1)
end subroutine out_of_range_error_limited_mixed_lowreal
!
subroutine out_of_range_error_limited_mixed_lowint
$(message,upper,ubound_in,lower,lbound_in)
implicit none
character(len=*), intent(in) :: message
double precision, intent(in) :: upper
integer , intent(in) :: lower
logical, intent(in) :: ubound_in
logical, intent(in) :: lbound_in
character(len=32) :: cupper
character(len=32) :: clower
character(len=1 ) :: c1u
character(len=1 ) :: c1l
write(out_channel,'(a)') ' Out of bounds error.'
write(out_channel,'(a)') message
write(out_channel,'(a)',advance='no') ' Legal range: '
if( lbound_in) c1l = '['
if(.not.lbound_in) c1l = '('
if( ubound_in) c1u = ']'
if(.not.ubound_in) c1u = ')'
write(cupper,*) upper
write(clower,*) lower
write(out_channel,*)
$ c1l,trim(adjustl(clower)),',',trim(adjustl(cupper)),c1u
call mrccend(1)
end subroutine out_of_range_error_limited_mixed_lowint
!
subroutine illegal_value_error_int
$( message , actual , expected )
character(len=*), intent(in) :: message
integer, intent(in) :: actual
integer, intent(in) :: expected
write(out_channel,'(a)') ' Illegal value error.'
write(out_channel,'(a)') message
write(out_channel,'(a)',advance='no') ' Expected:'
write(out_channel,*) expected
write(out_channel,'(a)',advance='no') ' Actual:'
write(out_channel,*) actual
call mrccend(1)
end subroutine illegal_value_error_int
!
subroutine illegal_value_error_char
$( message , actual , expected )
character(len=*), intent(in) :: message
character(len=*), intent(in) :: actual
character(len=*), intent(in) :: expected
write(out_channel,'(a)') ' Illegal value error.'
write(out_channel,'(a)') message
write(out_channel,'(a)',advance='no') ' Expected:'
write(out_channel,*) expected
write(out_channel,'(a)',advance='no') ' Actual:'
write(out_channel,*) actual
call mrccend(1)
end subroutine illegal_value_error_char
!
subroutine no_implementation_error( message )
character(len=*), intent(in) :: message
write(out_channel,'(a)') ' Expected feature is not implemented:'
write(out_channel,'(a)') message
write(out_channel,'(a)')
$' This state of the program will produce an undefined behavior.'
write(out_channel,'(2a)') ' The program will stop.'
call mrccend(1)
end subroutine no_implementation_error
!
end module error_handler
C
subroutine oslccversioncheck(localcc,scftype,mult,calc,iout)
implicit none
integer localccn,iout
character*4 localcc,mult
character*5 scftype
character*32 calc
logical lanycc
if (localcc.eq.'off ') then
localccn=0
else
read(localcc,*) localccn
endif
if(localccn.le.2021.and.localcc.ne.'off '.and.
$ ((scftype.ne.'rhf '.and.scftype.ne.' ').or.
$ (mult.ne.' '.and.mult.ne.'1 '))) then
if(localccn.lt.2022.and.lanycc(calc)) then
write(iout,"(a60)")
$ 'Error: Open-shell LNO-CC is only available via localcc>2021!'
call mrccend(1)
elseif(localccn.lt.2021.and..not.lanycc(calc)) then
write(iout,"(a60)")
$ 'Error: Open-shell LMP2 is only available via localcc>=2021!'
call mrccend(1)
endif
endif
end subroutine oslccversioncheck
Reference in New Issue View Git Blame Copy Permalink