API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Mar 20th, 2017
82
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 17.70 KB | None | 0 0
raw download clone embed print report
  1. MODULE fy_wavefuns
  2.  
  3.   USE types
  4.   USE cylho_basis_fixed
  5.  
  6.   IMPLICIT NONE
  7.  
  8.  
  9.   REAL(kind=r_kind), protected :: hbaromega0basis, hbaromegaperpbasis, hbaromegazbasis, eps2basis
  10.  
  11.   ! Transformation from asymptotic cyl. oscillator basis to canonical basis
  12.   ! |a>can = sum_alpha D(alpha,a) |alpha>cyl
  13.   REAL(kind=r_kind), protected, allocatable, target :: Dn(:,:), Dp(:,:)
  14.  
  15.   !Using rectangular matrices as number of canonical states < number of basis states.
  16.  
  17.  
  18. !!$  REAL(kind=r_kind), protected, target :: Dn(size_matels,size_matels), Dp(size_matels,size_matels)
  19.  
  20.   ! single-particle energies
  21. !!$  REAL(kind=r_kind), protected, target :: En(size_matels), Ep(size_matels)
  22.   REAL(kind=r_kind), protected, allocatable, target :: En(:), Ep(:)
  23.  
  24.   ! quasi-particle energies
  25.   REAL(kind=r_kind), protected, allocatable :: Eqpn(:), Eqpp(:)
  26.   INTEGER, protected, allocatable :: sort_qpn(:), sort_qpp(:)
  27.  
  28.   !BCS u and v factors in canonical basis
  29.   REAL(kind=r_kind), protected, allocatable :: ucann(:), vcann(:), ucanp(:), vcanp(:)
  30.  
  31.   INTEGER, protected :: no_neutrons, no_protons
  32.   REAL(kind=r_kind), protected :: efer_p, efer_n
  33.   REAL(kind=r_kind), protected :: efer2_p, efer2_n
  34.   REAL(kind=r_kind), protected :: Delta_p, G_p, Delta_n, G_n
  35.  
  36.  
  37.   !BCS U and V matrices in cylho basis
  38.   ! alpha(k)^dagger = ucan(k)a(k)^dagger -vcan(k)a(k_conj) = Sum_lambda u(k)D(lambda,k)a(lambda)^dagger - v(k)D(lambda,k_conj)^*a(lambda)
  39.   ! = Sum_lambda U(lambda,k)a(lambda)^dagger  + V(lambda,k)a(lambda)
  40.   COMPLEX(kind=r_kind), protected, allocatable :: Un(:,:), Vn(:,:), Up(:,:), Vp(:,:)
  41.  
  42.   COMPLEX(kind=r_kind), protected, allocatable :: rhon(:,:), kappan(:,:), rhop(:,:), kappap(:,:)
  43.  
  44.  
  45.  
  46.   INTEGER, protected :: no_statesn, no_statesp !number of canonical states not counting degeneracy = number of degenerate pairs
  47.   INTEGER, protected, allocatable, target :: sort_statesn(:), sort_statesp(:)
  48.   !INTEGER, protected, allocatable :: sort_parposstatesn(:), sort_parnegstatesn(:), sort_parposstatesp(:), sort_parnegstatesp(:)
  49.  
  50.  
  51.  
  52.   TYPE quasi_particle
  53.      INTEGER :: Omega2
  54.      INTEGER :: pi
  55.      INTEGER :: dom_comp !dominating component in unpairied canonical state
  56.      REAL(kind=r_kind) :: amp_dom_comp !amplitude of that comp
  57.      REAL(kind=r_kind) :: occ !v^2 for quasiparticle
  58.   END type quasi_particle
  59.  
  60.   type(quasi_particle), allocatable, protected, target :: quasiparticlesn(:), quasiparticlesp(:)
  61.  
  62.  
  63.   !TODO where to put phases for sigma = -1/2. Currently no phase is multiplied for sigma = -1/2 !done
  64.   !TODO block diagonality for canonical states !skip
  65.  
  66.   INTEGER, protected, target :: ompiblockposn(1:(2*nmax_matels+1),-1:1), ompiblocksizep(1:(2*nmax_matels+1),-1:1) ! omega2, pi
  67.   INTEGER, protected, target :: ompiblockposp(1:(2*nmax_matels+1),-1:1), ompiblocksizen(1:(2*nmax_matels+1),-1:1) ! omega2, pi
  68.  
  69.  
  70. CONTAINS
  71.  
  72.   !try to identify spherical shell of lowest energy quasiparticle
  73.   !by looking at nearly degenerate Eqps
  74.   SUBROUTINE print_spherical_shells(filename,filenameNordheim)
  75.     IMPLICIT NONE
  76.     CHARACTER(*), intent(in) :: filename, filenameNordheim
  77.     INTEGER :: II,JJ,dom_comp
  78.     INTEGER :: outunit, outunitNordheim
  79.     REAL(kind=r_kind) :: Eqp_min, Eqp
  80.     INTEGER :: Omeg2n_max, Omeg2p_max
  81.     INTEGER :: pin, pip, ln, lp, s2n, s2p, Jnp
  82.  
  83.     OPEN(file=TRIM(filename),newunit=outunit)
  84.  
  85.     WRITE(outunit,*) 'no_statesn = ',no_statesn
  86.     WRITE(outunit,*) 'neutron quasiparticles:'
  87.     WRITE(outunit,'(A4,A6,2A10,2A4,A6,2A10,4A4)') '#nr','qp_nr','Eqp','Esp','om2','pi','dc','amp','occ','N','nz','lam','om2'
  88.  
  89.     Eqp_min = Eqpn(2*sort_qpn(1)-1)
  90.     Omeg2n_max = ABS(quasiparticlesn(2*sort_qpn(1)-1)%omega2)
  91.     pin = quasiparticlesn(2*sort_qpn(1)-1)%pi
  92.  
  93.     DO II = 1,no_statesn
  94.  
  95.        JJ = sort_qpn(II)
  96.        dom_comp = quasiparticlesn(2*JJ-1)%dom_comp
  97.        
  98.        Eqp = Eqpn(2*JJ-1)
  99.        
  100.        IF( Eqp > Eqp_min + 0.1 .or. quasiparticlesn(2*JJ-1)%pi /= pin) CYCLE
  101.        
  102.        Omeg2n_max = MAX( Omeg2n_max, ABS( quasiparticlesn(2*JJ-1)%omega2 ) )
  103.  
  104.        WRITE(outunit,'(I4,I6,2F10.4,2I4,I6,2F10.4,4I4)') 2*II-1, 2*JJ-1, Eqpn(2*JJ-1),En(2*JJ-1), quasiparticlesn(2*JJ-1)%Omega2, quasiparticlesn(2*JJ-1)%pi, dom_comp , quasiparticlesn(2*JJ-1)%amp_dom_comp, quasiparticlesn(2*JJ-1)%occ, lookup_qn_matels(dom_comp)%n, lookup_qn_matels(dom_comp)%nz, lookup_qn_matels(dom_comp)%lambda, lookup_qn_matels(dom_comp)%omega2
  105.        dom_comp = quasiparticlesn(2*JJ)%dom_comp
  106.        WRITE(outunit,'(I4,I6,2F10.4,2I4,I6,2F10.4,4I4)') 2*II, 2*JJ, Eqpn(2*JJ),En(2*jj), quasiparticlesn(2*JJ)%Omega2, quasiparticlesn(2*JJ)%pi, dom_comp , quasiparticlesn(2*JJ)%amp_dom_comp, quasiparticlesn(2*JJ)%occ, lookup_qn_matels(dom_comp)%n, lookup_qn_matels(dom_comp)%nz, lookup_qn_matels(dom_comp)%lambda, lookup_qn_matels(dom_comp)%omega2
  107.  
  108.  
  109.     END DO
  110.  
  111.     WRITE(outunit,*)
  112.     WRITE(outunit,*) 'no_statesp = ',no_statesp
  113.     WRITE(outunit,*) 'proton quasiparticles:'
  114.     WRITE(outunit,'(A4,A6,2A10,2A4,A6,2A10,4A4)') '#nr','qp_nr','Eqp','Esp','om2','pi','dc','amp','occ','N','nz','lam','om2'
  115.    
  116.     Eqp_min = Eqpp(2*sort_qpp(1)-1)
  117.     Omeg2p_max = ABS(quasiparticlesp(2*sort_qpp(1)-1)%omega2)
  118.     pip = quasiparticlesp(2*sort_qpp(1)-1)%pi
  119.  
  120.     DO II = 1,no_statesp
  121.  
  122.        JJ = sort_qpp(II)
  123.        dom_comp = quasiparticlesp(2*JJ-1)%dom_comp
  124.  
  125.  
  126.        Eqp = Eqpp(2*JJ-1)
  127.        
  128.        IF( Eqp > Eqp_min + 0.1 .or. quasiparticlesp(2*JJ-1)%pi /= pip) CYCLE
  129.        
  130.        Omeg2p_max = MAX( Omeg2p_max, ABS( quasiparticlesp(2*JJ-1)%omega2 ) )
  131.  
  132.  
  133.        WRITE(outunit,'(I4,I6,2F10.4,2I4,I6,2F10.4,4I4)') 2*II-1, 2*JJ-1, Eqpp(2*JJ-1), Ep(2*JJ-1), quasiparticlesp(2*JJ-1)%Omega2, quasiparticlesp(2*JJ-1)%pi, dom_comp , quasiparticlesp(2*JJ-1)%amp_dom_comp, quasiparticlesp(2*JJ-1)%occ, lookup_qn_matels(dom_comp)%n, lookup_qn_matels(dom_comp)%nz, lookup_qn_matels(dom_comp)%lambda, lookup_qn_matels(dom_comp)%omega2
  134.        dom_comp = quasiparticlesp(2*JJ)%dom_comp
  135.        WRITE(outunit,'(I4,I6,2F10.4,2I4,I6,2F10.4,4I4)') 2*II, 2*JJ, Eqpp(2*JJ), Ep(2*JJ), quasiparticlesp(2*JJ)%Omega2, quasiparticlesp(2*JJ)%pi, dom_comp , quasiparticlesp(2*JJ)%amp_dom_comp, quasiparticlesp(2*JJ)%occ, lookup_qn_matels(dom_comp)%n, lookup_qn_matels(dom_comp)%nz, lookup_qn_matels(dom_comp)%lambda, lookup_qn_matels(dom_comp)%omega2
  136.  
  137.  
  138.     END DO
  139.  
  140.  
  141.  
  142.     CLOSE(outunit)
  143.  
  144.  
  145.     OPEN(file=TRIM(filenameNordheim),newunit=outunitNordheim)
  146.  
  147.     ln = (Omeg2n_max + 1)/2
  148.     s2n = 1
  149.     IF( (MOD(ln,2) == 0 .and. pin == 1) .or. (MOD(ln,2) == 1 .and. pin == -1)  ) THEN
  150.        !guessed correct ln
  151.     ELSE      
  152.        ln = (Omeg2n_max -1)/2
  153.        s2n = -1
  154.     END IF
  155.  
  156.     lp = (Omeg2p_max + 1)/2
  157.     s2p = 1
  158.     IF( (MOD(lp,2) == 0 .and. pip == 1) .or. (MOD(lp,2) == 1 .and. pip == -1)  ) THEN
  159.        !guessed correct lp
  160.     ELSE      
  161.        lp = (Omeg2p_max -1)/2
  162.        s2p = -1
  163.     END IF
  164.  
  165.     !Nordheim rule, intrisic spins parallel
  166.     IF( s2n == s2p) THEN
  167.        Jnp = (Omeg2n_max + Omeg2p_max)/2
  168.     ELSE
  169.        Jnp = ABS( (Omeg2n_max - Omeg2p_max)/2 )
  170.     END IF
  171.  
  172.  
  173.     WRITE(outunitNordheim,'(8A10)') '# ln','2*jn','pin','lp','2*jp','pip','Jnp','pinp'
  174.     WRITE(outunitNordheim,'(8I10)') ln, Omeg2n_max, pin, lp, Omeg2p_max, pip,Jnp,pin*pip
  175.  
  176.  
  177.     CLOSE(outunitNordheim)
  178.  
  179.  
  180.   END SUBROUTINE print_spherical_shells
  181.  
  182.  
  183.  
  184.   SUBROUTINE print_quasiparticles(filename)
  185.     IMPLICIT NONE
  186.  
  187.     CHARACTER(*), intent(in) :: filename
  188.     INTEGER :: II,JJ,dom_comp
  189.     INTEGER :: outunit
  190.  
  191.     OPEN(file=TRIM(filename),newunit=outunit)
  192.  
  193.     WRITE(outunit,*) 'no_statesn = ',no_statesn
  194.     WRITE(outunit,*) 'neutron quasiparticles:'
  195.     WRITE(outunit,'(A4,A6,2A10,2A4,A6,2A10,4A4)') '#nr','qp_nr','Eqp','Esp','om2','pi','dc','amp','occ','N','nz','lam','om2'
  196.     DO II = 1,no_statesn
  197.  
  198.        JJ = sort_qpn(II)
  199.        dom_comp = quasiparticlesn(2*JJ-1)%dom_comp
  200.        WRITE(outunit,'(I4,I6,2F10.4,2I4,I6,2F10.4,4I4)') 2*II-1, 2*JJ-1, Eqpn(2*JJ-1),En(2*JJ-1), quasiparticlesn(2*JJ-1)%Omega2, quasiparticlesn(2*JJ-1)%pi, dom_comp , quasiparticlesn(2*JJ-1)%amp_dom_comp, quasiparticlesn(2*JJ-1)%occ, lookup_qn_matels(dom_comp)%n, lookup_qn_matels(dom_comp)%nz, lookup_qn_matels(dom_comp)%lambda, lookup_qn_matels(dom_comp)%omega2
  201.        dom_comp = quasiparticlesn(2*JJ)%dom_comp
  202.        WRITE(outunit,'(I4,I6,2F10.4,2I4,I6,2F10.4,4I4)') 2*II, 2*JJ, Eqpn(2*JJ),En(2*jj), quasiparticlesn(2*JJ)%Omega2, quasiparticlesn(2*JJ)%pi, dom_comp , quasiparticlesn(2*JJ)%amp_dom_comp, quasiparticlesn(2*JJ)%occ, lookup_qn_matels(dom_comp)%n, lookup_qn_matels(dom_comp)%nz, lookup_qn_matels(dom_comp)%lambda, lookup_qn_matels(dom_comp)%omega2
  203.  
  204.  
  205.     END DO
  206.  
  207.     WRITE(outunit,*)
  208.     WRITE(outunit,*) 'no_statesp = ',no_statesp
  209.     WRITE(outunit,*) 'proton quasiparticles:'
  210.     WRITE(outunit,'(A4,A6,2A10,2A4,A6,2A10,4A4)') '#nr','qp_nr','Eqp','Esp','om2','pi','dc','amp','occ','N','nz','lam','om2'
  211.     DO II = 1,no_statesp
  212.  
  213.        JJ = sort_qpp(II)
  214.        dom_comp = quasiparticlesp(2*JJ-1)%dom_comp
  215.        WRITE(outunit,'(I4,I6,2F10.4,2I4,I6,2F10.4,4I4)') 2*II-1, 2*JJ-1, Eqpp(2*JJ-1), Ep(2*JJ-1), quasiparticlesp(2*JJ-1)%Omega2, quasiparticlesp(2*JJ-1)%pi, dom_comp , quasiparticlesp(2*JJ-1)%amp_dom_comp, quasiparticlesp(2*JJ-1)%occ, lookup_qn_matels(dom_comp)%n, lookup_qn_matels(dom_comp)%nz, lookup_qn_matels(dom_comp)%lambda, lookup_qn_matels(dom_comp)%omega2
  216.        dom_comp = quasiparticlesp(2*JJ)%dom_comp
  217.        WRITE(outunit,'(I4,I6,2F10.4,2I4,I6,2F10.4,4I4)') 2*II, 2*JJ, Eqpp(2*JJ), Ep(2*JJ), quasiparticlesp(2*JJ)%Omega2, quasiparticlesp(2*JJ)%pi, dom_comp , quasiparticlesp(2*JJ)%amp_dom_comp, quasiparticlesp(2*JJ)%occ, lookup_qn_matels(dom_comp)%n, lookup_qn_matels(dom_comp)%nz, lookup_qn_matels(dom_comp)%lambda, lookup_qn_matels(dom_comp)%omega2
  218.  
  219.  
  220.     END DO
  221.  
  222.  
  223.  
  224.     CLOSE(outunit)
  225.  
  226.   END SUBROUTINE print_quasiparticles
  227.  
  228.  
  229.  
  230.   SUBROUTINE setup_UV_matrices
  231.     IMPLICIT NONE
  232.  
  233.     INTEGER :: II
  234.     COMPLEX(kind=r_kind), allocatable :: testmat(:,:)
  235.     COMPLEX(kind=r_kind) :: tracerhon, tracerhop
  236.  
  237.     IF(ALLOCATED(Un)) DEALLOCATE(Un)
  238.     ALLOCATE(Un(size_matels,2*no_statesn))
  239.     IF(ALLOCATED(Vn)) DEALLOCATE(Vn)
  240.     ALLOCATE(Vn(size_matels,2*no_statesn))
  241.  
  242.     !these values should be overwritten
  243.     Un = -666.0
  244.     Vn = -666.0
  245.  
  246.     !note D is real
  247.     DO II = 1,2*no_statesn-1,2
  248.        Un(:,II) = Dn(:,II)*ucann(II)
  249.        Un(:,II+1) = Dn(:,II+1)*ucann(II)
  250.        Vn(:,II) = Dn(:,II+1)*(-1.0_r_kind)*vcann(II)
  251.        Vn(:,II+1) = Dn(:,II)*vcann(II)
  252.  
  253.        quasiparticlesn(II)%occ = vcann(II)**2
  254.        quasiparticlesn(II+1)%occ = vcann(II)**2
  255.  
  256.     END DO
  257.  
  258.  
  259.     IF(ALLOCATED(rhon)) DEALLOCATE(rhon)
  260.     ALLOCATE(rhon(size_matels,size_matels))
  261.     IF(ALLOCATED(kappan)) DEALLOCATE(kappan)
  262.     ALLOCATE(kappan(size_matels,size_matels))
  263.  
  264.     !rho = V^* V^T = V * V^T (V is real)
  265.     CALL ZGEMM('n','t',size_matels, size_matels, 2*no_statesn, one, Vn, size_matels ,Vn, size_matels, zero, rhon, size_matels)
  266.     !kappa = V^* U^T = V * U^T (V is real)
  267.     CALL ZGEMM('n','t',size_matels, size_matels, 2*no_statesn, one, Vn, size_matels ,Un, size_matels, zero, kappan, size_matels)
  268.  
  269.  
  270.     tracerhon = zero
  271.     DO II=1,size_matels
  272.        tracerhon = tracerhon + rhon(II,II)
  273.     END DO
  274.     WRITE(*,*) 'tracerhon = ',tracerhon
  275.  
  276.  
  277.     IF(ALLOCATED(Up)) DEALLOCATE(Up)
  278.     ALLOCATE(Up(size_matels,2*no_statesp))
  279.     IF(ALLOCATED(Vp)) DEALLOCATE(Vp)
  280.     ALLOCATE(Vp(size_matels,2*no_statesp))
  281.  
  282.     !these values should be overwritten
  283.     Up = -666.0
  284.     Vp = -666.0
  285.  
  286.     !note D is real
  287.     DO II = 1,2*no_statesp-1,2
  288.        Up(:,II) = Dp(:,II)*ucanp(II)
  289.        Up(:,II+1) = Dp(:,II+1)*ucanp(II)
  290.        Vp(:,II) = Dp(:,II+1)*(-1.0_r_kind)*vcanp(II)
  291.        Vp(:,II+1) = Dp(:,II)*vcanp(II)
  292.  
  293.        quasiparticlesp(II)%occ = vcanp(II)**2
  294.        quasiparticlesp(II+1)%occ = vcanp(II)**2
  295.  
  296.     END DO
  297.  
  298.     IF(ALLOCATED(rhop)) DEALLOCATE(rhop)
  299.     ALLOCATE(rhop(size_matels,size_matels))
  300.     IF(ALLOCATED(kappap)) DEALLOCATE(kappap)
  301.     ALLOCATE(kappap(size_matels,size_matels))
  302.  
  303.     !rho = V^* V^T = V * V^T (V is real)
  304.     CALL ZGEMM('n','t',size_matels, size_matels, 2*no_statesp, one, Vp, size_matels ,Vp, size_matels, zero, rhop, size_matels)
  305.     !kappa = V^* U^T = V * U^T (V is real)
  306.     CALL ZGEMM('n','t',size_matels, size_matels, 2*no_statesp, one, Vp, size_matels ,Up, size_matels, zero, kappap, size_matels)
  307.  
  308.  
  309.     tracerhop = zero
  310.     DO II=1,size_matels
  311.        tracerhop = tracerhop + rhop(II,II)
  312.     END DO
  313.     WRITE(*,*) 'tracerhop = ',tracerhop
  314.  
  315.  
  316.  
  317.  
  318. !!$    !Testing identities
  319. !!$    IF(ALLOCATED(testmat)) DEALLOCATE(testmat)
  320. !!$    ALLOCATE(testmat(2*no_statesn,2*no_statesn))
  321. !!$    testmat = -666.0
  322. !!$    testmat = MATMUL(TRANSPOSE(Dn),Dn)
  323. !!$
  324. !!$    WRITE(*,*) '(Dn^T*Dn)(1:10,1:10)'
  325. !!$    CALL zprint_matrix(10,testmat(1:10,1:10))
  326. !!$
  327. !!$    DO II = 1,2*no_statesn
  328. !!$       testmat(II,II) = testmat(II,II) - 1.0_r_kind
  329. !!$    END DO
  330. !!$    WRITE(*,*) 'Maxnorm of Dn^T * Dn -1 =',MAXVAL(ABS(testmat))
  331. !!$    testmat = -666.0
  332. !!$    testmat = MATMUL(TRANSPOSE(Un),Un) + MATMUL(TRANSPOSE(Vn),Vn)
  333. !!$    !subtract identity matrix
  334. !!$    DO II = 1,2*no_statesn
  335. !!$       testmat(II,II) = testmat(II,II) - 1.0_r_kind
  336. !!$    END DO
  337. !!$    WRITE(*,*) 'Maxnorm of Un^T * Un + Vn^T * Vn - 1 =',MAXVAL(ABS(testmat))
  338. !!$
  339. !!$    testmat = -666.0
  340. !!$    testmat = MATMUL(TRANSPOSE(Un),Vn) + MATMUL(TRANSPOSE(Vn),Un)
  341. !!$    WRITE(*,*) 'Maxnorm of Un^T * Vn + Vn^T * Un =',MAXVAL(ABS(testmat))
  342. !!$
  343. !!$    IF(ALLOCATED(testmat)) DEALLOCATE(testmat)
  344. !!$    ALLOCATE(testmat(size_matels,size_matels))
  345. !!$    testmat = -666.0
  346. !!$    testmat = MATMUL(Dn,TRANSPOSE(Dn))
  347. !!$    
  348. !!$    !FOR Dn*Dn^T, Un*Un^T, Vn*Vn^T:
  349. !!$    !ONLY THE ROWS AND COLUMNS WHICH CORRESPONDS TO BASIS STATES USED IN THE EXPNASION OF THE CANONICAL STATES WILL BE NONZERO!!
  350. !!$
  351. !!$    WRITE(*,*) '(Dn*Dn^T)(1:10,1:10)'
  352. !!$    CALL zprint_matrix(10,testmat(1:10,1:10))
  353. !!$    WRITE(*,*) '(Dn*Dn^T)(size-10:size,size-10:size)'
  354. !!$    CALL zprint_matrix(10,testmat(size_matels-10:size_matels,size_matels-10:size_matels))
  355. !!$
  356. !!$    DO II = 1,size_matels
  357. !!$       testmat(II,II) = testmat(II,II) - 1.0_r_kind
  358. !!$    END DO
  359. !!$    WRITE(*,*) 'Maxnorm of Dn * Dn^T -1 =',MAXVAL(ABS(testmat))
  360. !!$    testmat = -666.0
  361. !!$    testmat = MATMUL(Un,TRANSPOSE(Un)) + MATMUL(Vn,TRANSPOSE(Vn))
  362. !!$
  363. !!$    WRITE(*,*) '(Un*Un^T)(1:10,1:10)'
  364. !!$    CALL zprint_matrix(10,testmat(1:10,1:10))
  365. !!$    WRITE(*,*) '(Un*Un^T)(size-10:size,size-10:size)'
  366. !!$    CALL zprint_matrix(10,testmat(size_matels-10:size_matels,size_matels-10:size_matels))
  367. !!$
  368. !!$    !subtract identity matrix
  369. !!$    DO II = 1,size_matels
  370. !!$       testmat(II,II) = testmat(II,II) - 1.0_r_kind
  371. !!$    END DO
  372. !!$    WRITE(*,*) 'Maxnorm of Un * Un^T + Vn * Vn^T - 1 =',MAXVAL(ABS(testmat))
  373. !!$    testmat = -666.0
  374. !!$    testmat = MATMUL(Un,TRANSPOSE(Vn)) + MATMUL(Vn,TRANSPOSE(Un))
  375. !!$    WRITE(*,*) 'Maxnorm of Un * Vn^T + Vn * Un^T =',MAXVAL(ABS(testmat))
  376. !!$    IF(ALLOCATED(testmat)) DEALLOCATE(testmat)
  377.  
  378.  
  379.  
  380.   END SUBROUTINE setup_UV_matrices
  381.  
  382.  
  383.   SUBROUTINE print_matrix(size,matrix)
  384.  
  385.     USE types
  386.  
  387.     IMPLICIT NONE
  388.  
  389.     INTEGER, intent(in) :: size
  390.     REAL(kind=r_kind) :: matrix(size,size)
  391.  
  392.     INTEGER :: II, JJ
  393.  
  394.     DO II = 1,size
  395.        DO JJ = 1,size
  396.           WRITE(*,'(E18.8)',advance = "no") matrix(II,JJ)
  397.        END DO
  398.        WRITE(*,*)
  399.     END DO
  400.  
  401.  
  402.  
  403.   END SUBROUTINE PRINT_MATRIX
  404.  
  405.   SUBROUTINE zprint_matrix(size,matrix)
  406.  
  407.     USE types
  408.  
  409.     IMPLICIT NONE
  410.  
  411.     INTEGER, intent(in) :: size
  412.     COMPLEX(kind=r_kind) :: matrix(size,size)
  413.  
  414.     INTEGER :: II, JJ
  415.  
  416.     DO II = 1,size
  417.        DO JJ = 1,size
  418.           WRITE(*,'(A1,E9.2,A1,E9.2,A1)',advance = "no") "(",REAL(matrix(II,JJ),kind=r_kind),",",AIMAG(matrix(II,JJ)),")"
  419.        END DO
  420.        WRITE(*,*)
  421.     END DO
  422.  
  423.  
  424.  
  425.   END SUBROUTINE zprint_matrix
  426.  
  427.  
  428.  
  429.   SUBROUTINE canonical_uv_unpaired
  430.     IMPLICIT NONE
  431.  
  432.     INTEGER :: II
  433.     INTEGER :: last_occ_pair
  434.  
  435.     IF(ALLOCATED(ucann)) DEALLOCATE(ucann)
  436.     ALLOCATE(ucann(2*no_statesn))
  437.     IF(ALLOCATED(vcann)) DEALLOCATE(vcann)
  438.     ALLOCATE(vcann(2*no_statesn))
  439.  
  440.     IF(ALLOCATED(Eqpn)) DEALLOCATE(Eqpn)
  441.     ALLOCATE(Eqpn(2*no_statesn))
  442.  
  443.  
  444.     DO II = 1,2*no_statesn
  445.        IF(En(II) <= efer_n) THEN
  446.           ucann(II) = 0.0_r_kind
  447.           vcann(II) = 1.0_r_kind
  448.        ELSE
  449.           ucann(II) = 1.0_r_kind
  450.           vcann(II) = 0.0_r_kind
  451.        END IF
  452.  
  453.        Eqpn(II) = ABS(En(II)-efer_n)
  454.  
  455.     END DO
  456.  
  457.     IF(MOD(no_neutrons,2) == 1) THEN !equal filling approx for odd
  458.        last_occ_pair = no_neutrons/2
  459.        vcann(2*sort_statesn(last_occ_pair+1)-1) = 1/Sqrt(2.0_r_kind)
  460.        vcann(2*sort_statesn(last_occ_pair+1)) = 1/Sqrt(2.0_r_kind)
  461.        ucann(2*sort_statesn(last_occ_pair+1)-1) = 1/Sqrt(2.0_r_kind)
  462.        ucann(2*sort_statesn(last_occ_pair+1)) = 1/Sqrt(2.0_r_kind)
  463.  
  464.     END IF
  465.  
  466.     WRITE(*,*) 'Sum vcann**2 = ', DOT_PRODUCT(vcann,vcann)
  467.  
  468.     IF(ALLOCATED(ucanp)) DEALLOCATE(ucanp)
  469.     ALLOCATE(ucanp(2*no_statesp))
  470.     IF(ALLOCATED(vcanp)) DEALLOCATE(vcanp)
  471.     ALLOCATE(vcanp(2*no_statesp))
  472.  
  473.     IF(ALLOCATED(Eqpp)) DEALLOCATE(Eqpp)
  474.     ALLOCATE(Eqpp(2*no_statesp))
  475.  
  476.     DO II = 1,2*no_statesp
  477.        IF(Ep(II) <= efer_p) THEN
  478.           ucanp(II) = 0.0_r_kind
  479.           vcanp(II) = 1.0_r_kind
  480.        ELSE
  481.           ucanp(II) = 1.0_r_kind
  482.           vcanp(II) = 0.0_r_kind
  483.        END IF
  484.  
  485.        Eqpp(II) = ABS(Ep(II)-efer_p)
  486.  
  487.     END DO
  488.  
  489.     IF(MOD(no_protons,2) == 1) THEN !equal filling approx for odd
  490.        last_occ_pair = no_protons/2
  491.        vcanp(2*sort_statesp(last_occ_pair+1)-1) = 1/Sqrt(2.0_r_kind)
  492.        vcanp(2*sort_statesp(last_occ_pair+1)) = 1/Sqrt(2.0_r_kind)
  493.        ucanp(2*sort_statesp(last_occ_pair+1)-1) = 1/Sqrt(2.0_r_kind)
  494.        ucanp(2*sort_statesp(last_occ_pair+1)) = 1/Sqrt(2.0_r_kind)
  495.     END IF
  496.  
  497.     WRITE(*,*) 'Sum vcanp**2 = ', DOT_PRODUCT(vcanp,vcanp)
  498.  
  499.  
  500.  
  501.  
  502.   END SUBROUTINE canonical_uv_unpaired
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!