Skip to content
Hide whitespace changes
Inline Side-by-side
pfarm/cpu/src_mpi_omp_mol/energy_grid.f90 0 → 100644
View file @ f8829881
module energy_grid
! Define maximum scattering energy for energy grid region
! Time-stamp: "2003-09-19 13:19:02 cjn"
use precisn, only: wp
use io_units, only: fo
use mpi_params, only: io_processor
implicit none
logical, save :: fine_mesh ! fine/coarse energy grid
real(wp), save :: ebig ! largest scattering energy
private
public energy_max
public fine_mesh, ebig
contains
subroutine energy_max (imesh, nmesh)
! Define ebig, the maximum scattering energy for case
use rmx1_in, only: emax, targ_high_sce, uniform_e_grid, belthr
use hfile_data, only: ntarg, etarg, tchl, nc
integer, intent(in) :: imesh ! current mesh
integer, intent(in) :: nmesh ! number of energy regions
real(wp) :: e1
if (imesh == 1 .and. nmesh > 1) then ! resonance region (fine mesh)
fine_mesh = .true.
e1 = MINVAL(etarg(1:ntarg)) ! ground state energy threshold
if (targ_high_sce > 0) then ! sequence of final target
ebig = etarg(targ_high_sce) - e1 - belthr
else if (uniform_e_grid) then ! end on highest threshold
ebig = MAXVAL(etarg(1:ntarg)) - e1 - belthr
else ! end on highest theshold for this symmetry
ebig = MAXVAL(etarg(tchl(1:nc))) - e1 - belthr
end if
else ! coarse energy mesh
fine_mesh = .false.
ebig = emax
end if
if(io_processor) write (fo,'(/a,es12.4,a/)') 'Max. energy, ebig = ', ebig, &
' Scaled Rydbergs'
end subroutine energy_max
end module energy_grid
pfarm/cpu/src_mpi_omp_mol/error_out.f90 0 → 100644
View file @ f8829881
module error_out
implicit none
private
public error_check
contains
subroutine error_check (err, msg)
use io_units, only: fo
integer, intent(in) :: err ! flag
character(len=*), intent(in) :: msg ! error message
integer :: no
if (err /= 0) then
write (fo, fmt='(a)') msg
end if
end subroutine error_check
end module error_out
pfarm/cpu/src_mpi_omp_mol/io_units.f90 0 → 100644
View file @ f8829881
module io_units
! Time-stamp: "2003-06-12 07:37:43 cjn"
implicit none
private
public fi, fo, dsk, dskp1, dskp2
public set_dsk, set_dskp1, set_dskp2
! ang values:
public jbufiv, jbuffv, jbufie, jbuffe
public jbufiv1, jbuffv1, jbufiv2, jbuffv2
! ham values:
public jbuff1, jbuff2, jbuff4, jbuff5, itape1, itape3, jbufd
public jbufim, jbuffm, jbufev
public jibb, jfbb, jibc, jfbc, fl, ihunit, itm
public jbufivd, jbuffvd, jbufiv2d, jbuffv2d
public nfnml, fa
integer, parameter :: fi = 5 ! standard input
integer, parameter :: fo = 6 ! standard output
integer, parameter :: nfnml = 7 ! namelist file
integer, parameter :: fa = 9 ! sector data file
integer, save :: dsk = 0 ! scratch disk unit
integer, save :: dskp1 = 0 ! scratch disk unit
integer, save :: dskp2 = 0 ! scratch disk unit
integer, save :: jbufiv = 20
integer, save :: jbuffv = 21
integer, save :: jbuff1 = 24
integer, save :: jbuff2 = 25
integer, save :: itape1 = 26
integer, save :: jbuff4 = 27
integer, save :: jbufiv1 = 24 ! -> jbuff1
integer, save :: jbuffv1 = 25 ! -> jbuff2
integer, save :: jbufiv2 = 26 ! -> itape1
integer, save :: jbuffv2 = 27 ! -> jbuff4
integer, save :: jbuff5 = 28
integer, save :: jbufd = 29
integer, save :: itape3 = 10
integer, save :: jbufie = 22
integer, save :: jbuffe = 23
integer, save :: jbufim = 45
integer, save :: jbuffm = 46
integer, save :: jbufev = 47
! dipole units used in stmmat:
integer, save :: jbufivd = 41
integer, save :: jbuffvd = 42
integer, save :: jbufiv2d = 43
integer, save :: jbuffv2d = 44
integer, save :: jibb = 30 ! abbi
integer, save :: jfbb = 31 ! abbr
integer, save :: jibc = 32 ! abci
integer, save :: jfbc = 33 ! abcr
integer, save :: fl = 36 ! tgt dataset unit
integer, save :: ihunit = 37 ! unit for Hamiltonian output
integer, save :: itm = 38 ! unit for tarmom file
contains
subroutine set_dsk (dsk_unit)
integer, intent(in) :: dsk_unit
dsk = dsk_unit
return
end subroutine set_dsk
subroutine set_dskp1 (dsk_unit)
integer, intent(in) :: dsk_unit
dskp1 = dsk_unit
return
end subroutine set_dskp1
subroutine set_dskp2 (dsk_unit)
integer, intent(in) :: dsk_unit
dskp2 = dsk_unit
return
end subroutine set_dskp2
end module io_units
pfarm/cpu/src_mpi_omp_mol/mpi_params.f90 0 → 100644
View file @ f8829881
module mpi_params
implicit none
private
public taskid, numtasks
public io_processor
integer, save :: taskid, numtasks
logical, save :: io_processor=.false.
end module mpi_params
pfarm/cpu/src_mpi_omp_mol/pdg_ctl.f90 0 → 100644
View file @ f8829881
module pdg_ctl
! control diagonalization of sector Hamiltonians and RWA calculation
! Time-stamp: "2003-10-14 12:49:30 cjn"
! modified for jpi_coupling april to dec 2006 vmb
use precisn, only: wp
use io_units, only: fo
implicit none
private
public diag
contains
subroutine diag (nc, ncp, spins, sqno, imesh)
use hfile_data, only: lchl, kschl, tchl
use rw_amplitudes, only: start_sa_file, start_sa_file_id, close_sa_file, ampltd
use potl_cofs, only: cfs, cfsj, del_cf
use slp, only: qno
use def_sectors, only: nsect, nl, nx, w_sect, asect, reset_sec, &
sectors
use sector_potl, only: potl, potl_pams, reset_potl
use dist_mat, only: a_fill, sh_diag, kill_az
! use blacs, only: ctxt, p, q, myrow, mycol, io_processor, p0, &
! q0, p_error
use sec_ham, only: h_reset, gauleg, legndr, def_rb_vals, xi, def_ncq
use rmx1_in, only: inc_lrp_prop, packed_cf, bug9, jpi_coupling
use error_out, only: error_check
use mpi_params, only: io_processor, taskid, numtasks
! use mpi
integer, intent(in) :: nc ! # channels
integer, intent(in) :: ncp
type(qno), intent(in) :: sqno ! scattering q#s
integer, intent(in) :: spins(2) ! target spins
integer :: imesh ! COARSE,FINE mesh - imesh=1,2
real(wp), allocatable :: r(:)
real(wp) :: ral, rar
integer :: sect, st, status, nh, pp, qq
integer, save :: nsect_fine
integer :: nsect_coarse, mesh_taskid,coarse_mesh_first_taskid
integer ( kind = 4 ) ierr
call sectors (nc) ! define sector boundaries
! if (io_processor) then
! call start_sa_file (sqno, spins(1),sect)
! end if
call gauleg ! quadrature formula
call legndr ! Legendre functions
call potl_pams (nc, nx, ncp)
call def_ncq (ncp) ! pass ethr offset in h_el
! generate asymptotic potential coefficients
if (inc_lrp_prop) then
if (.not. packed_cf) then
if (.not. jpi_coupling) then
call cfs (nc, tchl(ncp:), lchl(ncp:), sqno, spins)
else
call cfsj (nc, tchl(ncp:), lchl(ncp:), kschl(ncp:), sqno, spins)
end if
end if
end if
allocate (r(nx), stat=status)
call error_check (status, 'pdiag: allocation r')
nh = nc * nl
write (fo,'(a,i0)') 'Total Number of Sectors = ', nsect
write (fo,'(a,i0)') 'Partition Hamiltonian dimension = ', nh
! Allow tasks to cycle through both meshes ensuring better load-balancing
! Ideally this would be achieved with a shared counter
if(imesh==1) then ! fine mesh - original taskid numbering for sectors
mesh_taskid = taskid
nsect_fine = nsect
else ! coarse mesh
nsect_coarse = nsect
mesh_taskid = taskid - nsect_fine
coarse_mesh_first_taskid = MOD(nsect_fine,numtasks)
mesh_taskid = MOD((taskid - coarse_mesh_first_taskid + numtasks),numtasks)
end if
do sect = 1, nsect
if (numtasks.gt.1) then
if (MOD(sect,numtasks)==MOD(mesh_taskid+1,numtasks)) then
if(imesh==1) then
write(fo,'(a,i0,a,i0,a,i0)')'FINE Region MPI Task ', taskid, ' calculating sector ', sect, ' , Hamiltonian dimension = ',nh
else
write(fo,'(a,i0,a,i0,a,i0)')'COARSE Region MPI Task ', taskid, ' calculating sector ', sect, ' , Hamiltonian dimension = ',nh
end if
call flush(fo)
else
cycle
end if
end if
call start_sa_file_id (sqno, spins(1), sect)
ral = asect(sect)
rar = asect(sect+1)
r = 0.5_wp * ((rar - ral) * xi + rar + ral)
call def_rb_vals (ral, rar)
call potl (nc, nx, r)
call A_fill (nh) ! form distributed sector H
call sh_diag (imesh, sect) ! diagonalize sector H
call ampltd (ral, rar, nc, nl, nh)
call kill_az
call close_sa_file()
end do
call h_reset
call reset_sec
call reset_potl
if ((inc_lrp_prop).and.(.not. packed_cf)) call del_cf
deallocate (r, stat=status)
call error_check (status, 'pdiag: deallocation')
end subroutine diag
end module pdg_ctl
pfarm/cpu/src_mpi_omp_mol/potl.f90 0 → 100644
View file @ f8829881
module sector_potl
! generate sector potential
! Time-stamp: "2003-09-19 11:52:40 cjn"
use precisn, only: wp
use io_units, only: fo
use potl_cofs, only: cf, lamd ! asymptotic potential coefficients
use rmx1_in, only: inc_lrp_prop, n_lambda, bug7, packed_cf
! use blacs, only: iam
implicit none
integer, save :: ion ! ionicity
integer, save :: lmx ! # multipoles
integer, allocatable, save :: lc(:) ! channel orbital a.m.
integer, save :: nc
real(wp), allocatable, save :: v(:,:,:)
private
public potl, potl_pams, reset_potl
public v, nc
contains
subroutine potl_pams (nct, nx, ncp)
! set potential parameters
use hfile_data, only: lchl
! use blacs, only: p_error
use error_out, only: error_check
use scaling, only: get_ionicity
integer, intent(in) :: nct ! # channels in partition
integer, intent(in) :: nx ! # radial points
integer, intent(in) :: ncp ! channel origin
integer :: status
nc = nct
allocate (v(nx,nc,nc), lc(nc), stat=status)
call error_check (status, 'potl_pams: allocation')
lc(1:nc) = lchl(ncp:ncp+nc-1)
ion = get_ionicity()
lmx = n_lambda
end subroutine potl_pams
subroutine reset_potl
! use blacs, only: ctxt, io_processor
use error_out, only: error_check
integer :: status
deallocate (v, lc, stat=status)
if (status /= 0) then
call error_check (status, 'reset_potl: allocation')
STOP
end if
end subroutine reset_potl
subroutine potl (nc, nx, r)
! use hfile_data, only: cf_nz, ic_label, ir_label, ncf_nonzero, ncf2_nonzero, spin_v, nc1
use hfile_data, only: cf_nz, ic_label, ir_label, ncf_nonzero, ncf2_nonzero, spin_v
! potl: potentials expanded in inverse powers of radial distance;
integer, intent(in) :: nc ! no. channels
integer, intent(in) :: nx ! no. basis fns
real(wp), intent(in) :: r(nx) ! radial points
real(wp) :: dg(nx)
integer :: ncf_temp(n_lambda), icf_st(n_lambda)
integer :: i, j, k, locc, lbb, lm, icf, kp1
if (inc_lrp_prop) then ! asymptotic potential included
lm = MIN(lmx, n_lambda)
v = 0.0_wp
if (packed_cf) then
! we just assume a symmetric matrix
if (spin_v == 1) then
icf_st(1:lm) = 1
ncf_temp(1:lm) = ncf_nonzero(1:lm)
else
icf_st(1:lm) = ncf_nonzero(1:lm) + 1
ncf_temp(1:lm) = ncf2_nonzero(1:lm)
end if
do j = 1, nc
dg = (REAL(lc(j)*(lc(j)+1),wp) / r - REAL(2*ion,wp)) / r
v(:,j,j) = dg
end do
do k = 1, lm
kp1 = k + 1
do icf = icf_st(k), ncf_temp(k)
j = ic_label(icf,k)
i = ir_label(icf,k)
v(:,i,j) = v(:,i,j) + cf_nz(icf,k) / r ** kp1
end do
end do
! symmetric matrix
do k = 1, lm
do icf = icf_st(k), ncf_temp(k)
j = ic_label(icf,k)
i = ir_label(icf,k)
v(:,j,i) = v(:,i,j)
end do
end do
else
! atomic method with symmetry taken into account
do j = 1, nc
dg = (REAL(lc(j)*(lc(j)+1),wp) / r - REAL(2*ion,wp)) / r
v(:,j,j) = dg
locc = j*(j-1)/2
do i = 1, j
do k = 1, lm, 2
lbb = (k + 1) / 2
! if (i == 1 .and. j == 3) &
! write(fo,'(a,2i3,2f14.7)') &
! '(1,3): lbb, lamd(lbb, rr(1), cf =', lbb, lamd(locc+i,lbb), &
! (1.0_wp/ r(1)), cf(locc+i,lbb)
v(:,i,j) = v(:,i,j) + cf(locc+i,lbb) &
/ r**(lamd(locc+i,lbb)+1)
v(:,j,i) = v(:,i,j)
! write(fo,*) 'iam', iam, ' i,j,k, lbb, cf:', i,j,k, lbb, cf(locc+i,lbb)
end do
end do
end do
end if
else
v = 0.0_wp
do j = 1, nc
dg = (REAL(lc(j)*(lc(j)+1),wp) / r - REAL(2*ion,wp)) / r
v(:,j,j) = dg
end do
end if
if (bug7 > 0) then
! write (fo,'(a)') ' radial values '
! write (fo,'(6f14.7)') r(:)
! do i = 1, nc
do i = 1, MIN(3,nc)
do j = 1, nc
! do j = 1, MIN(3,nc)
! write (fo,'(a,i4,i4)') ' Potential (v) matrix:, col, &
write (1025,'(a,i4,i4)') ' Potential (v) matrix:, col, &
&row ',i,j
write (1025,'(6e14.7)') v(:,i,j)
! write (fo,'(6f14.7)') v(:,i,j)
end do
end do
end if
end subroutine potl
end module sector_potl
pfarm/cpu/src_mpi_omp_mol/potl_coefs.f90 0 → 100644
View file @ f8829881
module potl_cofs
! Generate asymptotic coefficients from target multipoles
! Time-stamp: "2003-09-23 09:42:38 cjn"
! modified for jpi_coupling april to dec 2006 vmb
use precisn, only: wp
use io_units, only: fo
implicit none
private
public cfs, cfsj, cf, lamd, del_cf, rdtmom, rdjtmom
logical, save :: gams = .false.
integer, parameter :: ngam = 500 ! n entries in gam table
real(wp), save :: gam(ngam) ! gam table
type tmblock ! transition moment type definition
sequence
integer :: spn ! spin value
integer :: ntm ! # transition moments
integer, pointer :: lmds(:) ! multipole index
integer, pointer :: lbls(:) ! locator array
real(wp), pointer :: tms(:) ! transition moments
end type tmblock
integer, save :: ntmb ! # target spins
type(tmblock), allocatable, save :: tmb(:) ! target transition moments
real(wp), allocatable, save :: cf(:,:)
integer, allocatable, save :: lamd(:,:)
integer, save :: nt ! # targets
integer, allocatable, save :: tj(:), tl(:), tspn(:), tpty(:)
real(wp), allocatable, save :: tm(:,:)
contains
subroutine cfs (nc, tchl, lchl, sqno, spins)
! generate scaled asymptotic potential coefficients
use slp, only: qno, val_qno
use hfile_data, only: lamax, ltarg, starg, schl1, schl2, spin_v, cf_farm
! use blacs, only: ctxt, p_error
use scaling, only: scale_cf
use rmx1_in, only: farm_format
use error_out, only: error_check
integer, intent(in) :: nc ! no channels
integer, intent(in) :: tchl(nc) ! target sequence of channel
integer, intent(in) :: lchl(nc) ! orb. a.m. of channel
type(qno), intent(in) :: sqno ! scattering a. #s
integer, intent(in) :: spins(2) ! target spins required
real(wp) :: phz, z, drac
integer :: i, j, tari, tarj, li, lj, locc, tt, loct, tip
integer :: lbx, lbn, lb, lbb, tjl, tjs, tjp, til, tis
integer :: nspnt, qn(3), lrgl, nspn, npty, isp
integer :: status, smx, schl(nc), i_orig, j_orig
real(wp) :: a_val
smx = lamax
allocate (cf(nc*(nc+1)/2,(smx+1)/2), lamd(nc*(nc+1)/2,(smx+1)/2),&
stat=status)
call error_check (status, 'cfs: allocation error')
cf = 0.0_wp
lamd = 0
target_spins: do isp = 1, 2 ! 1 spin value or 2
nspnt = spins(isp)
if (nspnt == -999) cycle ! only one spin value
qn = val_qno(sqno)
lrgl = qn(1) ! scattering orb am
nspn = qn(2) ! scattering spin multiplicity
npty = qn(3) ! scattering parity
if (farm_format) then
! simplistic method: still copies several zeros
if (spin_v == 1) then
! spin_v updated by rmx1 for each call of pdiag to choose correct array
schl(1:nc) = schl1(1:nc)
else
schl(1:nc) = schl2(1:nc)
end if
do j = 1, nc
tarj = tchl(j)
! write(fo,'(a,i3,a,4i3)') 'iam', iam, ' j, tarj, starg(tarj), nspnt',j, tarj, starg(tarj), nspnt
if (starg(tarj) /= nspnt) cycle
lj = lchl(j)
tjl = ltarg(tarj)
locc = j*(j-1)/2
j_orig = schl(j)
do i = 1, j
tari = tchl(i)
! write(fo,'(a,i3,a,4i3)') 'iam', iam, ' i, tari, starg(tari), nspnt',i, tari, starg(tari), nspnt
if (starg(tari) /= nspnt) cycle
li= lchl(i)
til = ltarg(tari)
lbx = MIN(li+lj,til+tjl,smx) ! max lambda
lbn = MAX(ABS(li-lj),ABS(til-tjl),1) ! min lambda
i_orig = schl(i)
! if (i == j) write(fo,'(a,i3,a,4i3,2(1x,d14.6))') 'iam', iam, 'i,i_orig, j, j_orig, cf_farm, cf', &
! & i,i_orig, j, j_orig, cf_farm(i_orig,j_orig,1), cf(locc+i,1)
do lb = lbn, lbx, 1
lbb = (lb + 1) / 2
! write(fo,'(a,i3,a,4i3,2x,2i3,1x,d14.6)') 'iam', iam, 'i,i_orig, j, j_orig, lb, lbb, cf', &
! & i,i_orig, j, j_orig, lb, lbb, cf_farm(i_orig,j_orig,lb)
a_val = cf_farm(i_orig,j_orig,lb)
if (a_val /= 0.0_wp) then
cf(locc+i,lbb) = a_val
lamd(locc+i,lbb) = lb
end if
end do
end do
end do
cycle
end if
call get_tm (nspnt, smx, nt)
phz = 1
if (MOD(lrgl+npty,2) == 1) phz = -1
do j = 1, nc
tarj = tchl(j)
if (tspn(tarj) /= nspnt) cycle
lj = lchl(j)
tjl = tl(tarj)
tjs = tspn(tarj)
tjp = tpty(tarj)
locc = j*(j-1)/2
do i = 1, j
tari = tchl(i)
if (tspn(tari) /= nspnt) cycle
li= lchl(i)
til = tl(tari)
tis = tspn(tari)
tip = tpty(tari)
if (tari > tarj) then ! ensure input target sequencing
tt = tarj
tarj = tari
tari = tt
end if
loct = tarj * (tarj - 1) / 2 + tari
lbx = MIN(li+lj,til+tjl,smx) ! min lambda
lbn = MAX(ABS(li-lj),ABS(til-tjl),1) ! max lambda
do lb = lbn, lbx, 1
if (MOD(tjp+tip+lb,2) == 1) cycle
lbb = (lb + 1) / 2
if (ABS(tm(loct,lbb)) <= 1.0e-6_wp) cycle
drac = dracah(2*til, 2*tjl, 2*li, 2*lj, 2*lb, 2*lrgl)
z = phz * drac * rme(li, lj, lb)
cf(locc+i,lbb) = tm(loct,lbb) * z
lamd(locc+i,lbb) = lb
end do
end do
end do
deallocate (tm, stat=status)
if (status /= 0) STOP
end do target_spins
! cf_farm already scaled
if(.not. farm_format) call scale_cf (cf, lamd)
end subroutine cfs
subroutine cfsj (nc, tchl, lchl, kschl, sqno, spins)
! generate scaled asymptotic potential coefficients
use slp, only: qno, val_qno
use hfile_data, only: lamax, cf_farm
! use blacs, only: ctxt, p_error
use scaling, only: scale_cf
use rmx1_in, only: farm_format
use error_out, only: error_check
integer, intent(in) :: nc ! no channels
integer, intent(in) :: tchl(nc) ! target sequence of channel
integer, intent(in) :: lchl(nc) ! orb. a.m. of channel
integer, intent(in) :: kschl(nc) ! k q. no
type(qno), intent(in) :: sqno ! scattering q. #s
integer, intent(in) :: spins(2) ! k splits required
real(wp) :: phz, z, drac
integer :: i, j, tari, tarj, li, lj, locc, tt, loct, tip
integer :: lbx, lbn, lb, lbb, tjl, tjs, tjp, til, tis
integer :: ki, kj, kbar, tij, tjj
integer :: kspnt, qn(3), lrgl, nspn, npty, isp
integer :: status, smx
real(wp) :: a_val
smx = lamax
allocate (cf(nc*(nc+1)/2,(smx+1)/2), lamd(nc*(nc+1)/2,(smx+1)/2),&
stat=status)
call error_check (status, 'cfs: allocation error')
cf = 0.0_wp
lamd = 0
channel_k_splits: do isp = 1, 2 ! 1 k value or 2
kspnt = spins(isp) ! value of ksplit = 2k
if (kspnt == -999) cycle ! only one k value
qn = val_qno(sqno)
lrgl = qn(1) ! scattering orb am = 2J
nspn = qn(2) ! scattering spin multiplicity = 0
npty = qn(3) ! scattering parity
if (farm_format) then
! simple method: still copies several zeros
! nb if splitting is added, don't forget to add schl1/2 as in cfs
do j = 1, nc
tarj = tchl(j)
lj = lchl(j)
locc = j*(j-1)/2
do i = 1, j
tari = tchl(i)
li = lchl(i)
lbx = MIN(li+lj,smx) ! max lambda
lbn = MAX(ABS(li-lj),1) ! min lambda
do lb = lbn, lbx, 1
lbb = (lb + 1) / 2
a_val = cf_farm(i,j,lb)
if (a_val /= 0.0_wp) then
cf(locc+i,lbb) = a_val
lamd(locc+i,lbb) = lb
end if
end do
end do
end do
cycle
end if
call get_tm (0, smx, nt) ! target moments are not split
phz = 1
if (MOD(kspnt, 2) == 0) kbar=kspnt/2
if (MOD(kspnt, 2) == 1) kbar=(kspnt+1)/2! kbar remains real
if (MOD(npty+kbar, 2).eq.1) phz = -1
do j = 1, nc
kj = kschl(j)
if (kj /= kspnt) cycle
tarj = tchl(j)
lj = lchl(j)
tjj = tj(tarj)
tjl = tl(tarj)
tjs = tspn(tarj)
tjp = tpty(tarj)
locc = j*(j-1)/2
do i = 1, j
tari = tchl(i)
li = lchl(i)
ki = kschl(i)
if (ki /= kj) cycle
tij = tj(tari)
til = tl(tari)
tis = tspn(tari)
tip = tpty(tari)
if (tari > tarj) then ! ensure input target sequencing
tt = tarj
tarj = tari
tari = tt
end if
loct = tarj * (tarj - 1) / 2 + tari
lbx = MIN(li+lj,smx) ! max lambda
lbn = MAX(ABS(li-lj),1) ! min lambda
do lb = lbn, lbx, 1
if (MOD(tjp+tip+lb,2) == 1) cycle
lbb = (lb + 1) / 2
if (ABS(tm(loct,lbb)) <= 1.0e-6_wp) cycle
drac = dracah(tij, tjj, 2*li, 2*lj, 2*lb, kspnt)
z = phz * drac * rme(li, lj, lb)
cf(locc+i,lbb) = tm(loct,lbb) * z
lamd(locc+i,lbb) = lb
end do
end do
end do
deallocate (tm, stat=status)
if (status /= 0) STOP
end do channel_k_splits
! scaling already done for cf_farm
if(.not. farm_format) call scale_cf (cf, lamd)
end subroutine cfsj
subroutine del_cf
! use blacs, only: p_error
use error_out, only: error_check
integer :: status
deallocate (cf, lamd, stat=status)
call error_check (status, 'del_cf: deallocation error')
end subroutine del_cf
function dracah (i, j, k, l, m, n)
! racah coefficients
real(wp) :: dracah
integer, intent(in) :: i, j, k, l, m, n ! input momenta * 2
integer :: j1, j2, j3, j4, j5, j6, j7, numin, numax, icount, kk, ki
if (.not.gams) then
call gam_tab
gams = .true.
end if
j1 = i + j + m
j2= k + l + m
j3 = i + k + n
j4 = j + l + n
if (MOD(j1,2) /= 0 .or. MOD(j2,2) /= 0 .or. MOD(j3,2) /= 0 .or. &
MOD(j4,2) /= 0) then
dracah = 0.0_wp
return
end if
if (2*MAX(i,j,m) > j1 .or. 2*MAX(k,l,m) > j2 .or. 2*MAX(i,k,n) > j3 &
.or. 2*MAX(j,l,n) > j4) then
dracah = 0.0_wp
return
end if
j1 = j1 / 2
j2 = j2 / 2
j3 = j3 / 2
j4 = j4 / 2
j5 = (i + j + k + l) / 2
j6 = (i + l + m + n) / 2
j7 = (j + k + m + n) / 2
numin = MAX(j1, j2, j3, j4) + 1
numax = MIN(j5, j6, j7) + 1
dracah = 1.0_wp
icount = 0
if (numin < numax) then
do kk = numin+1, numax
ki = numax - icount
dracah = 1.0_wp - (dracah * real((ki * (j5-ki+2) * (j6-ki+2) * &
(j7-ki+2)),wp) / real((ki-1-j1) * (ki-1-j2) * (ki-1-j3) * &
(ki-1-j4),wp))
icount = icount + 1
end do
end if
dracah = dracah * (-1.0_wp)**(j5+numin+1) * &
EXP((gam(numin+1) - gam(numin-j1) - gam(numin-j2) - &
gam(numin-j3) - gam(numin-j4) - gam(j5+2-numin) - &
gam(j6+2-numin) - gam(j7+2-numin)) + ((gam(j1+1-i) + &
gam(j1+1-j) + gam(j1+1-m) - gam(j1+2) + gam(j2+1-k) + &
gam(j2+1-l) + gam(j2+1-m) - gam(j2+2) + gam(j3+1-i) + &
gam(j3+1-k) + gam(j3+1-n) - gam(j3+2) + gam(j4+1-j) + &
gam(j4+1-l) + gam(j4+1-n) - gam(j4+2)) / 2.0_wp))
end function dracah
function rme (l, lp, k)
! reduced matrix element (l//c(k)//lp)
! Ref: Fano and Racah, Irreducible Tensorial Sets, Ch 14, P. 81
real(wp) :: rme
integer, intent(in) :: l, lp ! initial, final orbital a.m.
integer, intent(in) :: k ! order of spherical tensor
integer :: i2g, ig, i1, i2, i3, imax, ier
real(wp) :: qusqrt
if (.not.gams) then
call gam_tab
gams = .true.
end if
if (k < ABS(l-lp) .or. k > l+lp) then
write (fo,'(a,/,a,2i3)') 'rme: a.m. &
& coupling error', 'rme zeroed for l, lp, k = ', l, lp, k
rme = 0.0_wp
return
end if
i2g = l + lp + k
ig = i2g / 2
if (MOD(i2g,2) /= 0) then
rme = 0.0_wp
return
end if
if (ig < 0) then
write (fo,'(a,/,a,2i3)') 'rme: &
& angle does not match', &
'rme zeroed for l, lp, k = ', l, lp, k
rme = 0.0_wp
return
elseif (ig == 0) then
rme = 1.0_wp
return
end if
i1 = ig - l
i2 = ig - lp
i3 = ig - k
imax = MAX(ig+1, i2g+2, 2*i3+1, 2*i2+1, 2*i1+1)
if (imax > ngam) then
write (fo,'(a,2i6)') 'rme: table overflow, imax, ngam = ', &
imax, ngam
STOP
end if
qusqrt = log(real(2*l+1,wp)) + log(real(2*lp+1,wp)) + gam(2*i1+1) +&
gam(2*i2+1) + gam(2*i3+1) - gam(i2g+2)
rme = EXP(0.5_wp * qusqrt + gam(ig+1) - gam(i1+1) - gam(i2+1) - &
gam(i3+1))
end function rme
subroutine gam_tab
! logs of factorials table
real(wp) :: x
integer :: i
gam(1) = 1.0_wp
gam(2) = 1.0_wp
x = 2.0_wp
do i = 3, 25
gam(i) = gam(i-1) * x
x = x + 1.0_wp
end do
gam(1:25) = LOG(gam(1:25))
x = 25.0_wp
do i = 26, ngam
gam(i) = gam(i-1) + LOG(x)
x = x + 1.0_wp
end do
end subroutine gam_tab
subroutine rdtmom
! read file of target moments
! use blacs, only: ictxt, p_error, io_processor
use io_units, only: fo, itm
use xdr_files, only: xdr_io, open_xdr, close_xdr
use rmx1_in, only: bug4, xdr_T_in, fltm
use error_out, only: error_check
real(wp), allocatable :: tmom(:,:)
integer, allocatable :: ntmom(:)
integer :: ibuf(2)
integer :: lamind, isp, ijl, ist, jst
integer :: lamlo, lamup, lam, lamn, nijl
integer :: ltmom, status, ios, s, i, j
integer :: nsv, spn, min_lam, max_lam, ntm
integer :: min_spn, max_spn, itmom, lamax, ip
if (xdr_T_in) then ! output file used XDR format
itmom = open_xdr (file=TRIM(fltm), action='read')
call xdr_io (itmom, ibuf, 2)
nt = ibuf(1)
lamax = ibuf(2)
allocate (tl(nt), tspn(nt), tpty(nt), stat=status)
if (status /= 0) then
write (fo,'(a,i6)') 'rdtmom: allocation error = ', status
STOP
end if
call xdr_io (itmom, tl, nt)
call xdr_io (itmom, tspn, nt)
call xdr_io (itmom, tpty, nt)
call xdr_io (itmom, ibuf, 2) ! upper, lower limits tgt spin
min_spn = ibuf(1)
max_spn = ibuf(2)
else ! fortran unformatted output file
itmom = itm
open (unit=itmom, file=TRIM(fltm), form='unformatted', &
status='old', action='read', iostat=ios)
if (ios /= 0) then
write (fo,'(a,i6)') 'rdtmom: error opening TARMOM, &
&iostat = ', ios
stop
end if
read (itmom) nt, lamax
allocate (tl(nt), tspn(nt), tpty(nt), stat=status)
call error_check (status, 'rdtmom: allocation error(1)')
read (itmom) tl(1:nt)
read (itmom) tspn(1:nt)
read (itmom) tpty(1:nt)
read (itmom) min_spn, max_spn
end if
lamind = (lamax + 1) / 2
ltmom = (nt*(nt+1)/2) * ((lamax+1)/2)
nsv = (max_spn - min_spn) / 2 + 1
allocate (tmom(ltmom,nsv), ntmom(nsv), stat=status)
if (status /= 0) then
write (fo,'(a,i6)') 'rdtmom: allocation error = ', status
STOP
end if
if (bug4 > 0) then ! echo data read from TARMOM
write (fo,'(a,i4,a,i4)') 'No. of states =', nt, &
' max lambda =', lamax
write (fo,'(a)') 'Target L:'
write (fo,'(25i3)') tl(1:nt)
write (fo,'(a)') 'Target S:'
write (fo,'(25i3)') tspn(1:nt)
write (fo,'(a)') 'Target P:'
write (fo,'(25i3)') tpty(1:nt)
write (fo,'(a)') 'Multiplicity limits:'
write (fo,'(25i3)') min_spn, max_spn
end if
s = 0
spins: do spn = min_spn, max_spn, 2
s = s + 1
if (xdr_T_in) then
call xdr_io (itmom, ibuf, 2)
isp = ibuf(1)
ntm = ibuf(2)
ntmom(s) = ntm
call xdr_io (itmom, tmom(1:ntm,s), ntm)
else
read (itmom) isp, ntm
ntmom(s) = ntm
read (itmom) tmom(1:ntm,s)
end if
if (bug4 > 0) then
write (fo,'(a,i4)') 'Target spin = ', isp
write (fo,'(a,i6)') 'Number of moments = ', ntm
write (fo,'(5f14.6)') tmom(1:ntm,s)
end if
end do spins
if (xdr_T_in) then
call close_xdr (itmom)
else
close (unit=itmom, status='keep')
end if
! unpack tmom:
ntmb = nsv
allocate (tmb(nsv), stat=status)
call error_check (status, 'rdtmom: tmb allocation error')
s = 0
do spn = min_spn, max_spn, 2
s = s + 1
ntm = ntmom(s)
allocate (tmb(s)%lbls(ntm), tmb(s)%lmds(ntm), tmb(s)%tms(ntm), &
stat=status)
call error_check (status, 'rdtmom: allocation error(5)')
tmb(s)%spn = spn
tmb(s)%ntm = ntm
ijl = 0
do i = 1, nt
if (tspn(i) /= spn) cycle
ip = i * (i - 1) / 2
do j = 1, i
if (tspn(j) /= spn) cycle
min_lam = ABS(tl(i) - tl(j))
if (MOD(tpty(i)+tpty(j)+min_lam,2) == 1) min_lam = &
min_lam + 1
if (min_lam == 0) min_lam = 2
max_lam = MIN(tl(i)+tl(j), lamax)
if ((min_lam+1)/2 <= (max_lam+1)/2) then
do lam = min_lam, max_lam, 2
ijl = ijl + 1
tmb(s)%lbls(ijl) = ip + j
tmb(s)%lmds(ijl) = lam
tmb(s)%tms(ijl) = tmom(ijl,s)
end do
end if
end do
end do
end do
deallocate (tmom, ntmom, stat=status)
call error_check (status, 'rdtmom: deallocation error')
end subroutine rdtmom
subroutine rdjtmom
! read file of target moments when jpi_coupling is true
! use blacs, only: ictxt, p_error, io_processor
use io_units, only: fo, itm
use xdr_files, only: xdr_io, open_xdr, close_xdr
use rmx1_in, only: bug4, xdr_T_in, fltm
use error_out, only: error_check
real(wp), allocatable :: tmom(:,:)
integer, allocatable :: ntmom(:)
integer :: ibuf(2)
integer :: lamind, isp, ijl, ist, jst
integer :: lamlo, lamup, lam, lamn, nijl
integer :: ltmom, status, ios, s, i, j
integer :: nsv, spn, min_lam, max_lam, ntm
integer :: min_spn, max_spn, itmom, lamax, ip
if (xdr_T_in) then ! output file used XDR format
itmom = open_xdr (file=TRIM(fltm), action='read')
call xdr_io (itmom, ibuf, 2)
nt = ibuf(1)
lamax = ibuf(2)
allocate (tj(nt), tl(nt), tspn(nt), tpty(nt), stat=status)
if (status /= 0) then
write (fo,'(a,i6)') 'rdjtmom: allocation error = ', status
STOP
end if
call xdr_io (itmom, tj, nt)
call xdr_io (itmom, tl, nt)
call xdr_io (itmom, tspn, nt)
call xdr_io (itmom, tpty, nt)
call xdr_io (itmom, ibuf, 2) ! upper, lower limits tgt spin
min_spn = ibuf(1)
max_spn = ibuf(2)
else ! fortran unformatted output file
itmom = itm
open (unit=itmom, file=TRIM(fltm), form='unformatted', &
status='old', action='read', iostat=ios)
if (ios /= 0) then
write (fo,'(a,i6)') 'rdjtmom: error opening TARMOM, &
&iostat = ', ios
stop
end if
read (itmom) nt, lamax
allocate (tj(nt), tl(nt), tspn(nt), tpty(nt), stat=status)
call error_check (status, 'rdjtmom: allocation error(1)')
read (itmom) tj(1:nt)
read (itmom) tl(1:nt)
read (itmom) tspn(1:nt)
read (itmom) tpty(1:nt)
read (itmom) min_spn, max_spn !set zero in fine.f90
end if
lamind = (lamax + 1) / 2
ltmom = (nt*(nt+1)/2) * ((lamax+1)/2)
nsv = (max_spn - min_spn) / 2 + 1 !nsv=1 here
allocate (tmom(ltmom,nsv), ntmom(nsv), stat=status)
if (status /= 0) then
write (fo,'(a,i6)') 'rdjtmom: allocation error = ', status
STOP
end if
if (bug4 > 0) then ! echo data read from TARMOM
write (fo,'(a,i4,a,i4)') 'No. of states =', nt, &
' max lambda =', lamax
write (fo,'(a)') 'Target 2J:'
write (fo,'(25i3)') tj(1:nt)
write (fo,'(a)') 'Target L:'
write (fo,'(25i3)') tl(1:nt)
write (fo,'(a)') 'Target S:'
write (fo,'(25i3)') tspn(1:nt)
write (fo,'(a)') 'Target P:'
write (fo,'(25i3)') tpty(1:nt)
write (fo,'(a)') 'Multiplicity limits:'
write (fo,'(25i3)') min_spn, max_spn
end if
s = 0
spins: do spn = min_spn, max_spn, 2! once thro loop
s = s + 1
if (xdr_T_in) then
call xdr_io (itmom, ibuf, 2)
isp = ibuf(1)
ntm = ibuf(2)
ntmom(s) = ntm
call xdr_io (itmom, tmom(1:ntm,s), ntm)
else
read (itmom) isp, ntm
ntmom(s) = ntm
read (itmom) tmom(1:ntm,s)
end if
if (bug4 > 0) then
write (fo,'(a,i4)') 'Notional Target spin = ', isp
write (fo,'(a,i6)') 'Number of moments = ', ntm
write (fo,'(5f14.6)') tmom(1:ntm,s)
end if
end do spins
if (xdr_T_in) then
call close_xdr (itmom)
else
close (unit=itmom, status='keep')
end if
! unpack tmom:
ntmb = nsv
allocate (tmb(nsv), stat=status)
call error_check (status, 'rdjtmom: tmb allocation error')
s = 0
do spn = min_spn, max_spn, 2
s = s + 1
ntm = ntmom(s)
allocate (tmb(s)%lbls(ntm), tmb(s)%lmds(ntm), tmb(s)%tms(ntm), &
stat=status)
call error_check (status, 'rdtmom: allocation error(5)')
tmb(s)%spn = spn
tmb(s)%ntm = ntm
ijl = 0
do i = 1, nt
ip = i * (i - 1) / 2
do j = 1, i
! min_lam = ABS(tl(i) - tl(j))
! if (MOD(tpty(i)+tpty(j)+min_lam,2) == 1) min_lam = &
! min_lam + 1
min_lam = MOD(tpty(i)+tpty(j),2)
if (min_lam == 0) min_lam = 2
! max_lam = MIN(tl(i)+tl(j), lamax)
max_lam = lamax
if ((min_lam+1)/2 <= (max_lam+1)/2) then
do lam = min_lam, max_lam, 2
ijl = ijl + 1
tmb(s)%lbls(ijl) = ip + j
tmb(s)%lmds(ijl) = lam
tmb(s)%tms(ijl) = tmom(ijl,s)
end do
end if
end do
end do
end do
deallocate (tmom, ntmom, stat=status)
call error_check (status, 'rdjtmom: deallocation error')
end subroutine rdjtmom
subroutine get_tm (spn, smx, nt)
! obtain expanded array of transition moments, tm for spin spn
use rmx1_in, only: bug5
integer, intent(in) :: spn ! target spin
integer, intent(in) :: smx ! max multipole
integer, intent(in) :: nt ! # targets
integer :: i, j, lp, ntm, l, ier, status, lb
do i = 1, ntmb
if (.not.associated(tmb(i)%tms)) exit
if (tmb(i)%spn /= spn) cycle ! not spin required
ntm = tmb(i)%ntm ! # transition moments for spin value
allocate (tm(nt*(nt+1)/2,(smx+1)/2), stat=status)
if (status /= 0) STOP
tm = 0.0_wp
do j = 1, ntm
l = tmb(i)%lmds(j)
lb = tmb(i)%lbls(j)
tm(lb, (l+1)/2) = tmb(i)%tms(j)
end do
if (bug5 > 0) then
do j = 1, ntm
l = tmb(i)%lmds(j)
lb = tmb(i)%lbls(j)
write (fo,'(3i6,a,f14.6)') j, lb, l, ' tms = ', tmb(i)%tms(j)
end do
end if
return ! required tm set has been located
end do
write (fo,'(a,i6)') 'get_tm: spin not found, spn = ', spn
STOP
end subroutine get_tm
end module potl_cofs
pfarm/cpu/src_mpi_omp_mol/precisn.f90 0 → 100644
View file @ f8829881
module precisn
! define precision-related paramters
! Time-stamp: "98/03/25 08:30:46 cjn"
implicit none
public
integer, parameter :: ibyte = 4
integer, parameter :: rbyte = 8
integer, parameter :: zbyte = 16
integer, parameter :: sp=selected_real_kind(6)
integer, parameter :: wp=selected_real_kind(12)
integer, parameter :: ep1=selected_real_kind(19)
integer, parameter :: ip_long = selected_int_kind(18)
integer, parameter :: ep = MAX(ep1, wp)
! real(wp), parameter :: acc8 = epsilon(1.0_wp)
! real(wp), parameter :: acc16 = epsilon(1.0_ep)
! real(wp), parameter :: fpmax = huge(1.0_wp) * acc8
! real(wp), parameter :: fpmin = tiny(1.0_wp) / acc8
! real(wp), save :: fplmax
! real(wp), save :: fplmin
real(wp), parameter :: acc8 = 2.0e-16_wp
real(wp), parameter :: acc16 = 3.0e-33_ep
real(wp), parameter :: fpmax = 1.0e60_wp
real(wp), parameter :: fpmin = 1.0e-60_wp
real(wp), parameter :: fplmax = 140.0_wp
real(wp), parameter :: fplmin = -140.0_wp
end module precisn