Commit 0ba38f36 by Matthias Redies

### remove d-notation

parent 45e3b727
 ... ... @@ -188,7 +188,7 @@ CONTAINS wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin) - usdus%us(l,n,jspin)*usdus%duds(l,n,jspin) IF (apw(l,n)) THEN fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin) dfj(l) = 0.0d0 dfj(l) = 0.0 ELSE dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk ... ...
 ... ... @@ -107,7 +107,7 @@ CONTAINS wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin)-usdus%us(l,n,jspin)*usdus%duds(l,n,jspin) !Wronski determinante IF (apw(l,n)) THEN fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin) dfj(l) = 0.0d0 dfj(l) = 0.0 ELSE dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk ... ...
 ... ... @@ -192,19 +192,19 @@ CONTAINS n_col = indxl2g(i, nb, hmat%blacsdata%mycol, 0, hmat%blacsdata%npcol) n_row = numroc (n_col, nb, hmat%blacsdata%myrow, 0, hmat%blacsdata%nprow) IF (hmat%l_real) THEN hmat%data_r(n_row+1:hmat%matsize1,i) = 0.d0 hmat%data_r(n_row+1:hmat%matsize1,i) = 0.0 ELSE hmat%data_c(n_row+1:hmat%matsize1,i) = 0.d0 hmat%data_c(n_row+1:hmat%matsize1,i) = 0.0 ENDIF ENDDO ! Use the ev_dist array to store the calculated values for the lower part. IF (hmat%l_real) THEN CALL pdtran(hmat%global_size1,hmat%global_size1,1.d0,hmat%data_r,1,1,& hmat%blacsdata%blacs_desc,0.d0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc) CALL pdtran(hmat%global_size1,hmat%global_size1,1.0,hmat%data_r,1,1,& hmat%blacsdata%blacs_desc,0.0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc) ELSE CALL pztranc(hmat%global_size1,hmat%global_size2,cmplx(1.d0,0.d0),hmat%data_c,1,1,& hmat%blacsdata%blacs_desc,cmplx(0.d0,0.d0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc) CALL pztranc(hmat%global_size1,hmat%global_size2,cmplx(1.0,0.0),hmat%data_c,1,1,& hmat%blacsdata%blacs_desc,cmplx(0.0,0.0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc) ENDIF ! Copy the calculated values to the lower part of the H matrix ... ... @@ -259,11 +259,11 @@ CONTAINS ! 2b. tmp2 = ev_dist**T IF (hmat%l_real) THEN CALL pdtran(ev_dist%global_size1,ev_dist%global_size1,1.d0,ev_dist%data_r,1,1,& ev_dist%blacsdata%blacs_desc,0.d0,tmp2_r,1,1,ev_dist%blacsdata%blacs_desc) CALL pdtran(ev_dist%global_size1,ev_dist%global_size1,1.0,ev_dist%data_r,1,1,& ev_dist%blacsdata%blacs_desc,0.0,tmp2_r,1,1,ev_dist%blacsdata%blacs_desc) ELSE CALL pztranc(ev_dist%global_size1,ev_dist%global_size1,cmplx(1.0,0.0),ev_dist%data_c,1,1,& ev_dist%blacsdata%blacs_desc,cmplx(0.d0,0.d0),tmp2_c,1,1,ev_dist%blacsdata%blacs_desc) ev_dist%blacsdata%blacs_desc,cmplx(0.0,0.0),tmp2_c,1,1,ev_dist%blacsdata%blacs_desc) ENDIF ! 2c. A = U**-T * tmp2 ( = U**-T * Aorig * U**-1 ) ... ... @@ -307,11 +307,11 @@ CONTAINS ! Set lower half from upper half IF (hmat%l_real) THEN CALL pdtran(hmat%global_size1,hmat%global_size1,1.d0,hmat%data_r,1,1,& hmat%blacsdata%blacs_desc,0.d0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc) CALL pdtran(hmat%global_size1,hmat%global_size1,1.0,hmat%data_r,1,1,& hmat%blacsdata%blacs_desc,0.0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc) ELSE CALL pztranc(hmat%global_size1,hmat%global_size1,cmplx(1.0,0.0),hmat%data_c,1,1,& hmat%blacsdata%blacs_desc,cmplx(0.d0,0.d0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc) hmat%blacsdata%blacs_desc,cmplx(0.0,0.0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc) ENDIF ... ... @@ -396,11 +396,11 @@ CONTAINS ! mult_ah_b_complex needs the transpose of U**-1, thus tmp2 = (U**-1)**T IF (hmat%l_real) THEN CALL pdtran(smat%global_size1,smat%global_size1,1.d0,smat%data_r,1,1,& smat%blacsdata%blacs_desc,0.d0,tmp2_r,1,1,smat%blacsdata%blacs_desc) CALL pdtran(smat%global_size1,smat%global_size1,1.0,smat%data_r,1,1,& smat%blacsdata%blacs_desc,0.0,tmp2_r,1,1,smat%blacsdata%blacs_desc) ELSE CALL pztranc(smat%global_size1,smat%global_size1,cmplx(1.d0,0.d0),smat%data_c,1,1,& smat%blacsdata%blacs_desc,cmplx(0.d0,0.d0),tmp2_c,1,1,smat%blacsdata%blacs_desc) CALL pztranc(smat%global_size1,smat%global_size1,cmplx(1.0,0.0),smat%data_c,1,1,& smat%blacsdata%blacs_desc,cmplx(0.0,0.0),tmp2_c,1,1,smat%blacsdata%blacs_desc) ENDIF #if defined (CPP_ELPA_201705003) ... ...
 ... ... @@ -119,8 +119,8 @@ ENDDO IF ( l_mcd ) THEN ! create an energy grid for mcd-spectra e_lo = 9.9d+9 e_up = -9.9d+9 e_lo = 9.9*10.0**9 e_up = -9.9*10.0**9 DO jspin = 1,input%jspins DO n = 1,atoms%ntype DO icore = 1 , mcd%ncore(n) ... ...
 ... ... @@ -12,11 +12,11 @@ module m_corespec ! PARAMETERS complex, parameter :: cone = cmplx(1.d0,0.d0) complex, parameter :: cimu = cmplx(0.d0,1.d0) real, parameter :: alpha = 7.29735257d-3 real, parameter :: mec2 = 0.51099891d6 real, parameter :: ecoredeep = 0.5d0 complex, parameter :: cone = cmplx(1.0,0.0) complex, parameter :: cimu = cmplx(0.0,1.0) real, parameter :: alpha = 7.29735257*10.0**-3 real, parameter :: mec2 = 0.51099891*10.0**6 real, parameter :: ecoredeep = 0.5 integer, parameter :: edgel(11) = (/0,1,1,2,2,3,3,4,4,5,5/) integer, parameter :: edgej(11) = (/1,1,3,3,5,5,7,7,9,9,11/) ... ...
 ... ... @@ -42,10 +42,10 @@ MODULE m_corespec_io csi%edge = "" csi%edgeidx = 0 csi%lx = -1 csi%ek0 = 0.d0 csi%emn = -2.d0 csi%emx = 20.d0 csi%ein = 0.1d0 csi%ek0 = 0.0 csi%emn = -2.0 csi%emx = 20.0 csi%ein = 0.1 csi%nqphi = 12 csi%nqr = 20 ... ... @@ -182,12 +182,12 @@ MODULE m_corespec_io &"maximum l: ","csi%lx = ",csi%lx,"will be used" ! csi%ek0 if(csi%ek0.le.0.d0) then if(csi%ek0.le.0.0) then write(*,csmsgs) trim(smeno),"found csi%ek0 <= 0.0 !"//csmsgerr ; stop endif csi%ek0 = csi%ek0*1000.d0 ! conversion from keV to eV csv%gamma = 1.d0+csi%ek0/mec2 csv%beta = sqrt(1.d0-1.d0/(csv%gamma**2)) csi%ek0 = csi%ek0*1000.0 ! conversion from keV to eV csv%gamma = 1.0+csi%ek0/mec2 csv%beta = sqrt(1.0-1.0/(csv%gamma**2)) if(csi%verb.eq.1) then write(*,csmsgses) trim(smeno),& &"kinetic energy of incoming electrons: ","csi%ek0 = ",csi%ek0,& ... ... @@ -204,8 +204,8 @@ MODULE m_corespec_io if(csi%emn.gt.csi%emx) then write(*,csmsgs) trim(smeno),"found csi%emn > csi%emx !"//csmsgerr ; stop endif if(csi%ein.le.0.d0) then write(*,csmsgs) trim(smeno),"found csi%ein <= 0.d0 !"//csmsgerr ; stop if(csi%ein.le.0.0) then write(*,csmsgs) trim(smeno),"found csi%ein <= 0.0 !"//csmsgerr ; stop endif if(((csi%emx-csi%emn)/csi%ein)-int((csi%emx-csi%emn)/csi%ein).ne.0) then write(*,csmsgs) trim(smeno),& ... ... @@ -216,7 +216,7 @@ MODULE m_corespec_io csv%egrid = (/(csi%emn+csi%ein*dble(i), i = 0,csv%nex)/) csv%nen = 0 !!\$ do i = 0,csv%nex !!\$ if(csv%egrid(i).ge.0.d0) then !!\$ if(csv%egrid(i).ge.0.0) then !!\$ csv%nen = i !!\$ exit !!\$ endif ... ... @@ -240,9 +240,9 @@ MODULE m_corespec_io &csv%nen,"will be used" if(.not.allocated(csv%eedge)) allocate(csv%eedge(csv%nljc)) csv%eedge = 0.d0 csv%eedge = 0.0 if(.not.allocated(csv%occ)) allocate(csv%occ(csv%nljc)) csv%occ = 0.d0 csv%occ = 0.0 l_cs = .true. ... ...
 ... ... @@ -44,7 +44,7 @@ CONTAINS gg = rk(k)*gb(l) IF ( apw(l) ) THEN fj(k,l,ispin) = 1.0*con1 * ff / us(l,ispin) gj(k,l,ispin) = 0.0d0 gj(k,l,ispin) = 0.0 ELSE IF (l_flag) THEN DO jspin = 1, jspins ... ... @@ -167,7 +167,7 @@ CONTAINS gg = lapw%rk(k,intspin)*gb(l) IF ( apw(l) ) THEN fj(k,l,ispin,intspin) = 1.0*con1 * ff / usdus%us(l,n,ispin) gj(k,l,ispin,intspin) = 0.0d0 gj(k,l,ispin,intspin) = 0.0 ELSE IF (noco%l_constr.or.l_socfirst) THEN DO jspin = 1, input%jspins ... ...
 ... ... @@ -398,8 +398,8 @@ SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_s !---> update overlap and l-diagonal hamiltonian matrix kj_end = MIN(ki,lapw%nv(iintsp)) VecHelpS = 0.d0 VecHelpH = 0.d0 VecHelpS = 0.0 VecHelpH = 0.0 DO l = 0,atoms%lmax(n) fjkiln = fj(ki,l,jintsp) gjkiln = gj(ki,l,jintsp) ... ...
 ... ... @@ -88,7 +88,7 @@ CONTAINS gg = lapw%rk(k,iintsp)*gb(l) ! IF ( apw(l) ) THEN ! fj(k,l,n,iintsp) = 1.0*con1 * ff / usdus%us(l,n,isp) ! gj(k,l,n,iintsp) = 0.0d0 ! gj(k,l,n,iintsp) = 0.0 ! ELSE !---> in a spin-spiral calculation fj and gj are needed !---> both interstitial spin directions at the same time ... ...
 ... ... @@ -139,7 +139,7 @@ CONTAINS wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin) - usdus%us(l,n,jspin)*usdus%duds(l,n,jspin) IF (apw(l,n)) THEN fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin) dfj(l) = 0.0d0 dfj(l) = 0.0 ELSE dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk ... ...
 ... ... @@ -297,7 +297,7 @@ else IF (i1.EQ.1) nn = 0 IF (i1.EQ.2) nn = nsz(1) zhelp2(:,:) = 0.d0 zhelp2(:,:) = 0.0 DO j = 1,nsize DO i = 1,nsz(jsp) zhelp2(i,j) = CONJG(hso(i+nn,j)) ... ... @@ -305,11 +305,11 @@ else ENDDO ! j if (l_real) THEN CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd,CMPLX(1.d0,0.d0),CMPLX(zmat(jsp)%data_r(:,:)),& zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(1,1,jsp2),zmat(1)%matsize1) CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd,CMPLX(1.0,0.0),CMPLX(zmat(jsp)%data_r(:,:)),& zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(1,1,jsp2),zmat(1)%matsize1) else CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd, CMPLX(1.d0,0.d0),zmat(jsp)%data_c(:,:),& zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(:,:,jsp2),zmat(1)%matsize1) CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd, CMPLX(1.0,0.0),zmat(jsp)%data_c(:,:),& zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(:,:,jsp2),zmat(1)%matsize1) endif ENDDO !isp ... ...
 ... ... @@ -287,10 +287,10 @@ c------------------------------------------------------------------- w=0 nqvect=0 nshort=0 ReJq=0.d0 ImJq=0.d0 ReJq=0.0 ImJq=0.0 sqsin=(sin(thetaJ))**2 tpi = 2.d0 * pimach() tpi = 2.0 * pimach() limit=nmagn-1 IF (nmagn.gt.mtypes) limit=mtypes ... ... @@ -375,7 +375,7 @@ c... DO ii=1,2 Dabsq(:)=ABS(q(:,nn,qcount)+((-1)**ii)*q(:,nnn,qcount)) IDabsq(:)=NINT(Dabsq(:)) divi(:)=ABS(Dabsq(:)/FLOAT(IDabsq(:))-1.d0) divi(:)=ABS(Dabsq(:)/FLOAT(IDabsq(:))-1.0) IF(((Dabsq(1).LT.tol).OR.(divi(1).LT.tol)).AND. & ((Dabsq(2).LT.tol).OR.(divi(2).LT.tol)).AND. & ((Dabsq(3).LT.tol).OR.(divi(3).LT.tol)))THEN ... ... @@ -397,12 +397,12 @@ c... c... Now calculate Jq=Re(Jq)+i*Im(Jq) mu=1 DO imt=1,mtypes ReJq(mu,mu,qcount)=-2.d0*(seigv(mu,mu,qcount,1) ReJq(mu,mu,qcount)=-2.0*(seigv(mu,mu,qcount,1) & -seigv0(mu,mu,1))/(M(mu)*M(mu)*sqsin) ImJq(mu,mu,qcount)=0.d0 ImJq(mu,mu,qcount)=0.0 DO remt=mu+1,mu+nmagtype(imt)-1 ReJq(remt,remt,qcount)=ReJq(mu,mu,qcount) ImJq(remt,remt,qcount)=0.d0 ImJq(remt,remt,qcount)=0.0 ENDDO!remt mu=mu+nmagtype(imt) ENDDO !imt ... ... @@ -411,10 +411,10 @@ c... Now calculate Jq=Re(Jq)+i*Im(Jq) DO nu=mu+1,nmagn ReJq(mu,nu,qcount)=((seigv0(mu,nu,2)- & seigv(mu,nu,qcount,1))/(M(mu)*M(nu)*sqsin)) & -(0.5d0*M(mu)*ReJq(mu,mu,qcount)/M(nu))- & (0.5d0*M(nu)*ReJq(nu,nu,qcount)/M(mu)) & -(0.5*M(mu)*ReJq(mu,mu,qcount)/M(nu))- & (0.5*M(nu)*ReJq(nu,nu,qcount)/M(mu)) IF(invs)THEN ImJq(mu,nu,qcount)=0.d0 ImJq(mu,nu,qcount)=0.0 ELSE ImJq(mu,nu,qcount)=((seigv(mu,nu,qcount,2) & -seigv(mu,nu,qcount,1))/ ... ... @@ -479,7 +479,7 @@ c ... for one magnetic atom per unit cell qcount=nqpt-1 lwork=2*nshort ALLOCATE (Cmat(qcount,nshort),DelE(qcount),work(lwork)) Cmat=0.d0 Cmat=0.0 IF (nshort.GE.nqpt)THEN WRITE(*,*) ' Please supply the data for', nshort, & 'q-points different from zero' ... ... @@ -492,7 +492,7 @@ c ... for one magnetic atom per unit cell scp=(q(1,1,n)*R(1,atsh,nn) & +q(2,1,n)*R(2,atsh,nn) & +q(3,1,n)*R(3,atsh,nn))*tpi Cmat(n,nn)=Cmat(n,nn)-1.d0+cos(scp) Cmat(n,nn)=Cmat(n,nn)-1.0+cos(scp) ENDDO ENDDO DelE(n)=ReJq(1,1,n)*2000 ! multiply by 2000 to get [mRy/muB**2] ... ... @@ -504,7 +504,7 @@ c ... for one magnetic atom per unit cell & work,lwork,info) c The routine dgels returns the solution, J(n), in the array DelE(n) Tc=0.d0 Tc=0.0 DO n=1,nshort Tc=Tc+nat(n)*DelE(n) !Mean-field Tc=1/3*(Sum_i(J_0,i)) WRITE(115,5005) n,lenR(n),DelE(n) ! J in units [mRy/muB**2] ... ... @@ -527,7 +527,7 @@ c... Perform the back-Fourier transform wrJ=0 DO atsh=1,nat(nnn) IF(atsh.gt.shmax) STOP 'jcoff2:increase shmax!' J=0.d0 J=0.0 DO n=1,nqpt-1 DO nn=1,nop IF(w(nn,n).EQ.1)THEN ... ... @@ -539,7 +539,7 @@ c... Perform the back-Fourier transform ENDIF ENDDO !nn ENDDO !n (qpts) J=(J/float(nqvect))*2000.d0 ! J in units [mRy/muB**2] J=(J/float(nqvect))*2000.0 ! J in units [mRy/muB**2] DO i=1,wrJ !A check for non-equivalent sub-shells IF(ABS(J-Jw(i)).LE.(tol))GOTO 55 ENDDO ... ... @@ -565,7 +565,7 @@ c... In case of only one magnetic atom per unit cell, calculate the mean-field mu=mu+nmagtype(imt) ENDDO !imt IF(nmagn.EQ.1) THEN Tc=157.889*M(1)*M(1)*Tc/3.d0 Tc=157.889*M(1)*M(1)*Tc/3.0 WRITE(115,*) '# Tc(mean field)= ',Tc ENDIF 5008 FORMAT(i4,i4,7(1x,f14.10)) ... ...
 ... ... @@ -44,7 +44,7 @@ SUBROUTINE gen_map(atoms,sym,oneD,hybrid) ratom = 0 DO ieq1 = 1,atoms%neq(itype) IF( all(abs(modulo(rtaual-atoms%taual(:,iatom0 + ieq1)+1d-12,1d0)).lt. 1d-10) ) THEN IF( all(abs(modulo(rtaual-atoms%taual(:,iatom0 + ieq1)+10.0**-12,1.0)).lt. 10.0**-10) ) THEN ratom = iatom0 + ieq1 hybrid%map ( iatom,isym) = ratom hybrid%tvec(:,iatom,isym) = nint ( rtaual-atoms%taual(:,ratom) ) ... ...
 ... ... @@ -107,7 +107,7 @@ C WRITE (ibfile,'('' sum = dvec**2 = '',f13.6)') sum ENDDO sum = sqrt(sum) ddist(n) = 0.5d0*sum ddist(n) = 0.5*sum WRITE (ibfile,'(/'' ddist('',i3,'')=(.5*sum**.5) '',f13.6)') > n,ddist(n) sum = 1.0/sum ... ...
 ... ... @@ -82,8 +82,8 @@ C C ---> save and data statements c save one,zero,half,eps,eps1 data zero/0.0d0/,one/1.0d0/,half/0.5d0/, + eps/1.0d-8/,eps1/1.0d-5/ data zero/0.0/,one/1.0/,half/0.5/, + eps/1.0*10.0**-8/,eps1/1.0*10.0**-5/ c c----------------------------------------------------------------------- if (kpri .ge. 3) then ... ...
 ... ... @@ -96,7 +96,7 @@ IF( all( matmul(rot(:,:,i),rot(:,:,j)) & .eq.reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))) & .and.all(modulo(matmul(rot(:,:,i),rtau(:,j))+rtau(:,i),1.0) & .lt.1d-10) )THEN & .lt.10.0**-10) )THEN IF(invtab(i).ne.0) STOP 'kptgen: inverse operation & & already defined.' invtab(i) = j ... ... @@ -137,7 +137,7 @@ DO k = 2,nsym rarr = matmul(rrot(:,:,k),bk(:,i)) * kpts%nkpt3 iarr2 = nint(rarr) IF(any(abs(iarr2-rarr).gt.1d-10)) THEN IF(any(abs(iarr2-rarr).gt.10.0**-10)) THEN WRITE(6,'(A,I3,A)') 'kptgen: Symmetry operation',k, & ' incompatible with k-point set.' ldum = .true. ... ... @@ -280,7 +280,7 @@ c function modulo1 maps kpoint into first BZ nkpt3 = (/a,b,c/) modulo1 = kpoint*nkpt3 help = nint(modulo1) IF(any(abs(help-modulo1).gt.1d-8)) THEN IF(any(abs(help-modulo1).gt.10.0**-8)) THEN modulo1 = kpoint*nkpt3 WRITE(*,*) modulo1 help = nint(modulo1) ... ... @@ -291,7 +291,7 @@ c function modulo1 maps kpoint into first BZ + 'modulo1: argument not an element of k-point set.', + calledby = 'kptgen_hybrid:modulo1') END IF modulo1 = modulo(help,nkpt3)*1d0/nkpt3 modulo1 = modulo(help,nkpt3)*1.0/nkpt3 END FUNCTION modulo1 ... ...
 ... ... @@ -138,9 +138,9 @@ C do 2200 i=1,3 dist=dist+(kvc(i,j)-kvc(i,l))**2 2200 continue if ( dist.lt.(dm*(1.0d0-1.0d-6)) ) dm=dist if ( dist.lt.(dm*(1.0-10.0**-6)) ) dm=dist 2100 continue if ( d.lt.(dm*(1.0d0+1.0d-6)) ) goto 2000 if ( d.lt.(dm*(1.0+10.0**-6)) ) goto 2000 d=dm lmin=l 2000 continue ... ... @@ -164,9 +164,9 @@ C do 3300 i=1,3 dist=dist+(kvc(i,j)-knew(i))**2 3300 continue if ( dist.lt.(dm*(1.0d0-1.0d-6)) ) dm=dist if ( dist.lt.(dm*(1.0-10.0**-6)) ) dm=dist 3200 continue if ( dm.lt.(dmax*(1.0d0+1.0d-6)) ) goto 3000 if ( dm.lt.(dmax*(1.0+10.0**-6)) ) goto 3000 dmax=dm do 3400 i=1,3 kc(i)=knew(i) ... ... @@ -177,7 +177,7 @@ C WE HAVE FOUND A NEW POINT KC . IF THE MINIMAL DISTANCE OF C THIS POINT IS LESS THAN THE DISTANCE D WE ALREADY HAVE , C WE STOP THE DYNAMIC CYCLE . C if ( dmax.lt.1.0001d0*d ) goto 4000 if ( dmax.lt.1.0001*d ) goto 4000 do 3500 i=1,3 kvc(i,lmin)=kc(i) 3500 continue ... ... @@ -188,7 +188,7 @@ C WE HAVE FOUND A REASONABLE DISTRIBUTION . TO AVOID TETRAHEDRA C WITH A FUNNY SHAPE , WE PUT ALL THE POINTS WHICH ARE NEAR A C SIDE ON THIS SIDE . C d=0.49d0*sqrt(d) d=0.49*sqrt(d) do 4100 l=ncorn+1,nkpt xmin=1.0e6 do n1=1,ncorn-2 ... ...
 ... ... @@ -108,8 +108,8 @@ C C ---> save and data statements c save one,zero,half,eps,eps1,iplus,iminus data zero/0.0d0/,one/1.0d0/,half/0.5d0/, + eps/1.0d-8/,eps1/1.0d-5/,iplus/1/,iminus/-1/ data zero/0.0/,one/1.0/,half/0.5/, + eps/1.0*10.0**-8/,eps1/1.0*10.0**-5/,iplus/1/,iminus/-1/ c c----------------------------------------------------------------------- c ... ...
 ... ... @@ -43,9 +43,9 @@ C Fred Hutson C C************************************************************************* C data eps/1.0d-10/ data eps1/1.0d-5/ pi = 4.0d0*atan(1.0d0) data eps/10.0**-10/ data eps1/10.0**-5/ pi = 4.0*atan(1.0) C C CONSTRUCT THE FIRST TETRAHEDRON C ... ... @@ -53,10 +53,10 @@ C C C FIND K-POINT NEAREST TO 1 C dm=1.0d9 dm=10.0**9 i2=0 do 100 i=2,nkpt dist=0.0d0 dist=0.0 do 110 icom=1,3 dist=dist+(kvc(icom,1)-kvc(icom,i))**2 110 continue ... ... @@ -71,19 +71,19 @@ C C FIND POINT NEAREST TO (1+I2)/2 , NOT ON THE LINE C CONNECTING 1 AND I2 C dm=1.0d9 dm=10.0**9 i3=0 do 200 i=2,nkpt if ( i.eq.i2 ) goto 200 dl=0.0d0 dist=0.0d0 dl=0.0 dist=0.0 do 210 icom=1,3 vect=kvc(icom,i)-0.5d0*(kvc(icom,1)+kvc(icom,i2)) vect=kvc(icom,i)-0.5*(kvc(icom,1)+kvc(icom,i2)) dist=dist+vect*vect dl=dl+vect*(kvc(icom,1)-kvc(icom,i2)) 210 continue dl=dl/dnorm if ( abs(dist-dl*dl).lt.0.01d0*dist ) goto 200 if ( abs(dist-dl*dl).lt.0.01*dist ) goto 200 if ( dist.gt.dm ) goto 200 dm=dist i3=i ... ... @@ -100,31 +100,31 @@ C \$ (kvc(1,1)-kvc(1,i2))*(kvc(3,i2)-kvc(3,i3)) cn(3)=(kvc(1,1)-kvc(1,i2))*(kvc(2,i2)-kvc(2,i3))- \$ (kvc(2,1)-kvc(2,i2))*(kvc(1,i2)-kvc(1,i3)) dnorm=0.0d0 dnorm=0.0 do 300 icom=1,3 dnorm=dnorm+cn(icom)**2 300 continue dnorm=1.0d0/sqrt(dnorm) dnorm=1.0/sqrt(dnorm) do 310 icom=1,3 cn(icom)=cn(icom)*dnorm 310 continue i4=0 dm=1.0d9 dm=10.0**9 do 400 i=2,nkpt if ( (i.eq.i2).or.(i.eq.i3) ) goto 400 dist=0.0d0 dl=0.0d0 dist=0.0 dl=0.0 do 410 icom=1,3 vect=kvc(icom,i)-(kvc(icom,1)+kvc(icom,i2)+ \$ kvc(icom,i3))/3.0d0 \$ kvc(icom,i3))/3.0 dist=dist+vect**2 dl=dl+vect*cn(icom) 410 continue if ( dl*dl.lt.0.01d0*dist ) goto 400 if ( dl*dl.lt.0.01*dist ) goto 400 if ( dist.gt.dm ) goto 400 i4=i dm=dist vt=dl/(dnorm*6.0d0) vt=dl/(dnorm*6.0) 400 continue IF ( i4==0 ) CALL juDFT_error(" tetcon3 ",calledby ="tetcon") ntetra(4,1)=i4 ... ... @@ -185,12 +185,12 @@ C \$ (kcorn(1,2)-kcorn(1,1))*(kcorn(3,3)-kcorn(3,1)) norm(3)=(kcorn(1,2)-kcorn(1,1))*(kcorn(2,3)-kcorn(2,1))- \$ (kcorn(2,2)-kcorn(2,1))*(kcorn(1,3)-kcorn(1,1)) vol=0.0d0 vol=0.0 do 1400 i=1,3 vol=vol+norm(i)*(kvc(i,ntetra(j,it))-kcorn(i,1)) 1400 continue vol=vol/6.0d0 if ( vol.lt.0.0d0 ) goto 1500 vol=vol/6.0 if ( vol.lt.0.0 ) goto 1500 do 1600 i=1,3 norm(i)=-norm(i) 1600 continue ... ... @@ -201,10 +201,10 @@ C STORE THE K-POINT ADDRESS IN NKADD ARRAY ACCORDING TO THE C ORDER OF DISTANCE BETWEEN THE K-POINT AND THIS FACE. C do 1800 nkp=1,nkpt length=0.0d0 length=0.0 do 1850 i=1,3 length=length+(kvc(i,nkp) & -(kcorn(i,1)+kcorn(i,2)+kcorn(i,3))/3.0d0)**2 & -(kcorn(i,1)+kcorn(i,2)+kcorn(i,3))/3.0)**2 1850 continue ddist(nkp)=length nkadd(nkp)=nkp ... ... @@ -235,11 +235,11 @@ C kcorn(i,4)=kvc(i,nkadd(nkp)) 2050 continue nside(4)=nkadd(nkp) vt=0.0d0 vt=0.0 do 2100 i=1,3 vt=vt+norm(i)*(kcorn(i,4)-kcorn(i,1)) 2100 continue vt=vt/6.0d0 vt=vt/6.0 C C REJECT POINT NKP IF IT IS ON THE WRONG SIDE . C ... ... @@ -253,7 +253,7 @@ C C FIRST CHECK IF TETRAHEDRON NK2 IS ON THE DANGEROUS SIDE C do 3150 i=1,4 sum=0.0d0 sum=0.0 do 3160 jj=1,3 sum=sum+norm(jj)*(kvc(jj,ntetra(i,nk2))-kcorn(jj,1)) 3160 continue ... ... @@ -268,7 +268,7 @@ C do i=1,4 do 3200 jj=1,4 if(ntetra(i,nk2).eq.nside(jj)) goto 3200 sum=0.0d0 sum=0.0 do icom=1,3 cn(icom)=kvc(icom,ntetra(i,nk2))-kcorn(icom,jj) sum=sum+cn(icom)*cn(icom) ... ... @@ -277,11 +277,11 @@ C cn(icom)=cn(icom)/sqrt(sum) end do do k=1,l sum=0.0d0 sum=0.0 do icom=1,3 sum=sum+cn(icom)*d(icom,k) end do if((1.0d0-sum).lt.eps) goto 3200 if((1.0-sum).lt.eps) goto 3200 C C EXCLUDE THE VECTOR WHICH HAS THE SAME DIRECTION AS AN EXISTING ONE. C ... ... @@ -312,19 +312,19 @@ C & .and.(abs(dnm(3)).lt.eps)) goto 3500 do 3600 k=1,l if((k.eq.i).or.(k.eq.jj)) goto 3600 sum=0.0d0 sum=0.0 do icom=1,3 sum=sum+d(icom,k)*dnm(icom) end do if(abs(sum).lt.eps) goto 3600 do 3700 kk=1,l if((kk.eq.i).or.(kk.eq.jj).or.(kk.eq.k)) goto 3700 sss=0.0d0 sss=0.0 do icom=1,3 sss=sss+d(icom,kk)*dnm(icom) end do