!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Fortran90 programs to verify the formulas given in the paper ! ! ! ! "Analytic Quadratic Integration Over the ! ! Two-Dimensional Brillouin Zone" ! ! ! ! Frank E. Harris ! ! ! ! This program last revised 10 November 2001 ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! These programs are not intended for production purposes, but for ! validation and demonstration of the method. ! The actual integration formulas are in the subroutine entitled ! "tri_int", which includes the logic for identifying the integration ! paths as well as evaluating the integrals. This subroutine is ! reasonably efficient and is probably suitable as a part of a ! production program. ! ! The subroutine "triangles" subdivides an initial triangle to the ! extent indicated in its input, creating an array identifying the ! integration points needed for each triangle. ! ! The subroutines "v2wts" and "e2q" convert between function values ! at individual points and coefficients in an expansion through ! quadratic terms. ! ! The subroutine "weights" uses the subprograms described above to ! generate a set of weights based on input energy (and fermi energy) ! data ! ! No attempt has been made to ensure numerical stability on extremely ! fine grids ! ! The subprograms "fcalc", "ecalc", and "load_pts" are of value only ! for test purposes: Respectively, they (1) provide functional forms ! for functions to be integrated, (2) provide energy profiles from ! which the integration region is specified, and (3) load these data ! for use in test calculations. ! ! The main program, "integral_test", takes as input a fermi energy, ! and then runs a large number of tests at a series of subdivisions ! of the regions of integration, using the subroutine "test_int". ! The output displays the area within the fermi surface (which is ! exact for some energy profiles, approximate for others) and the ! value of the quantity being integrated (exact if quadratic or less ! and over an exactly represented profile, approximate otherwise). ! Tests are set up for the irreducible wedge of a BZ with P4 symmetry ! and for the entire BZ. ! ! The test cases set up here are for five different functions to be ! integrated: (1) r**2 (will integrate exactly if profile is exact), ! (2) r (will integrate approximately), (3) cos(2 pi x) + cos(2 pi y), ! (4) cos(6 pi x) + cos(6 pi y), and (5) 0.2 x**3 + 0.3 y**3. ! ! Three energy profiles have been set up: (1) E = r**2 (makes circular ! arcs which can be represented exactly), (2) r**2/(1+0.2 cos(4 phi), ! where phi is the angle in a polar coordinate system (this makes a ! contour that must be approximated), and (3) 0.4(x-0.65)**4 +0.2(y-0.33)**4, ! which will produce a region of integration that must be approximated ! and is located away from the origin. The integration region is where ! E < E_fermi. ! ! -------------------------------------------------------------------------- program integral_test implicit none real*8 :: E integer*4 :: k write(6,*) 'Input Fermi energy for test:' read(5,*) E do k=1,30 call test_int(k,E) enddo stop end ! -------------------------------------------------------------------------- subroutine test_int(k,E) ! Makes table of integral values at grids of subdivision 2**n (n=0..5) in ! each dimension for test problem k at Fermi energy E ! If k <= 15 (wedge), uses a single triangular grid with vertices at ! (x,y) = (0,0), (1,0), and (1,1). If k > 15 (BZ), combines results from ! two grids: first for triangle with vertices (-1,-1), (1,-1), and ! (1,1); then for triangle with vertices (1,1), (-1,1), and (-1,-1). implicit none integer*4 :: n,i,j,npts,ntri,k real*8 :: e,pi,w,f,phi,area character*23 :: text(2) integer*4,allocatable :: ipoint(:,:) real*8,allocatable :: wts(:),ee(:),ff(:) data text / 'E < e_fermi in wedge ', 'E < e-fermi in whole BZ'/ pi = 4.d0*atan(1.d0) write(6,991) k,E 991 format('Test no. ',i2,' Fermi energy =',f12.8) i = mod(k-1,5)+1 j = mod((k-1)/5,3)+1 n = (k-1)/15+1 write(6,997) i,j,text(n) 997 format('F no.',i2,' E profile no.',i2,3X,A23) write(6,*) ' N Area F integral' do n=0,5 npts = (2**(n+1)+1)*(2**n+1) ntri = 2**(2*n) allocate(ipoint(6,ntri),wts(npts),ee(npts),ff(npts)) call triangles(ipoint,n) call load_pts(ff,ee,n,k) call weights(ee,e,n,ipoint,wts) w = 0.d0 f = 0.d0 do i=1,npts w = w + wts(i) f = f + wts(i)*ff(i) ! write(6,902) i,wts(i) 902 format(' Point no. and weight ',i3,f12.8) enddo if (k > 15) then call load_pts(ff,ee,n,-k) call weights(ee,e,n,ipoint,wts) do i = 1,npts w = w + wts(i) f = f + wts(i)*ff(i) enddo w = 4.d0*w ! scale adjustment f = 4.d0*f endif write(6,903) n,w,f 903 format(i6,2f14.8) deallocate(ipoint,wts,ee,ff) enddo return end ! -------------------------------------------------------------------------- function ecalc(x,y,j) ! Energy profiles 1, 2, and 3 real*8 :: ecalc,x,y,phi integer*4 :: j go to (1,2,3),j 1 ecalc = x**2 + y**2 return 2 phi = atan2(y,x) ecalc = (x**2+y**2)/(1.d0+0.2d0*cos(4.d0*phi)) return 3 ecalc = 5.d0*(x-0.65d0)**4+2.d0*(y-0.33d0)**4 return end ! -------------------------------------------------------------------------- function fcalc(x,y,j) ! Functions to be integrated 1 through 5 real*8 :: fcalc,x,y,pi,s integer*4 :: j pi = 4.d0*atan(1.d0) go to (1,2,3,4,5),j 1 fcalc = x**2 + y**2 return 2 fcalc = sqrt(x**2 + y**2) return 3 s = 2.d0*pi fcalc = cos(s*x) + cos(s*y) return 4 s = 6.d0*pi fcalc = cos(s*x) + cos(s*y) return 5 fcalc = 0.2d0*x**3 + 0.3d0*y**3 return end ! -------------------------------------------------------------------------- subroutine load_pts(ff,ee,n,k) ! Loads arrays ff and ee with values for all triangle-index ! points where the x (and y) intervals are subdivided into ! subintervals of length 1/2**(n+1) the original interval, for ! x >= y. Identifies points using same scheme as "ipt" (used ! in "triangles"), obtains ff and ee values from external ! functions fcalc and ecalc respectively. ! Assumes the first 15 test cases to be for the irreducible ! wedge of a square BZ, a triangle with vertices at (0,0), ! (1,0), and (1,1). The remaining test cases are for one ! half of the entire BZ: if k>0 the triangle with vertices (x,y)= ! (-1,-1), (1,-1), and (1,1). If k<0 the triangle with vertices ! (1,1), (-1,1), and (-1,-1). (These two triangles are combined ! in the test calculation.) implicit none integer*4 :: n,it,ix,iy,k,ij(2,30),kk real*8 :: ff(*),ee(*),fcalc,ecalc,delta,x,y,xy0 data ij / 1,1, 2,1, 3,1, 4,1, 5,1, 1,2, 2,2, 3,2, 4,2, 5,2, & 1,3, 2,3, 3,3, 4,3, 5,3, 1,1, 2,1, 3,1, 4,1, 5,1, & 1,2, 2,2, 3,2, 4,2, 5,2, 1,3, 2,3, 3,3, 4,3, 5,3 / kk = abs(k) if (abs(k) .le. 15) then delta = 1.d0/2**(n+1) xy0 = 0.d0 else if (k .gt. 15) then delta = 1.d0/2**n xy0 = -1.d0 else delta = -1.d0/2**n xy0 = 1.d0 endif it = 0 do ix = 0,2**(n+1) x = xy0 + ix*delta do iy = 0,ix it = it + 1 y = xy0 + iy*delta ff(it) = fcalc(x,y,ij(1,kk)) ee(it) = ecalc(x,y,ij(2,kk)) enddo enddo end ! -------------------------------------------------------------------------- subroutine weights(ee,e,n,ipoint,wts) ! For previously loaded set of ee (energy profile) and input e ! value (the fermi energy), using input triangle data, obtain ! combined weights (summed over all triangles) for quadratic ! integration, with wts(it) containing weight for point it ! (using triangular indexing scheme). Grid subdivision 2**n ! in each dimension. Output is scaled to grid of total area 0.5. implicit none integer*4 :: n,it,ipoint(6,*) real*8 :: ee(*),e,wts(*),e_tri(6),q(6),v(6),wt_tri(6) do it=1,(2**(n+1)+1)*(2**n+1) wts(it) = 0.d0 enddo do it = 1,2**(2*n) e_tri = ee(ipoint(:,it)) call e2q(e_tri,q) call tri_int(e,q,v,1) ! option 1: normal operation call v2wts(v,wt_tri) wts(ipoint(:,it)) = wts(ipoint(:,it)) + wt_tri enddo do it=1,(2**(n+1)+1)*(2**n+1) wts(it) = wts(it)/2**(2*n) enddo end ! -------------------------------------------------------------------------- subroutine v2wts(v,wt_tri) ! Converts set of 6 monomial integrals for the occupied ! portion of a triangle to a set of weights for second- ! order integration for the occupied region using the ! six points of the triangle implicit none integer*4 :: i,j,MM(6,6) real*8 :: v(6),wt_tri(6) data MM/ 1,-3,-3,2,4,2, 0,-1,0,2,0,0, 0,0,-1,0,0,2, & 0,4,0,-4,-4,0, 0,0,0,0,4,0, 0,0,4,0,-4,-4 / do j = 1,6 wt_tri(j) = 0.d0 do i = 1,6 wt_tri(j) = wt_tri(j) + v(i)*MM(i,j) enddo enddo end ! -------------------------------------------------------------------------- subroutine e2q(e_tri,q) ! Converts set of 6 e_tri values for a triangle to ! coefficients q for corresponding quad expansion implicit none integer*4 :: i,j,MM(6,6) real*8 :: e_tri(6),q(6) data MM/ 1,-3,-3,2,4,2, 0,-1,0,2,0,0, 0,0,-1,0,0,2, & 0,4,0,-4,-4,0, 0,0,0,0,4,0, 0,0,4,0,-4,-4 / do i = 1,6 q(i) = 0.d0 do j = 1,6 q(i) = q(i) + MM(i,j)*e_tri(j) enddo enddo end ! -------------------------------------------------------------------------- subroutine triangles(ipoint,n) ! Sets up point numbering for points 1 thru 6 of ! each triangle when each side of the original triangle ! is subdivided into 2**n subtriangle sides [making the ! total number of triangles 2**(2n)]. Array ipoint(pt_no, ! tri_no) contains the points of triangle tri_no using the ! triangular indexing scheme in function ipt. The triangle ! numbering starts from 1, the six points of each triangle ! are in the order (for the standard simplex) (0,0),(1,0), ! (0,1),(1/2,0),(1/2,1/2),(0,1/2). implicit none integer*4 :: n,ipoint(6,*),it,ix,iy it = 0 do ix = 1,2**n do iy = 0,ix if (iy .ne. ix) then it = it + 1 ipoint(1,it) = ipt(2*ix,2*iy) ipoint(2,it) = ipt(2*ix,2*iy+2) ipoint(3,it) = ipt(2*ix-2,2*iy) ipoint(4,it) = ipt(2*ix,2*iy+1) ipoint(5,it) = ipt(2*ix-1,2*iy+1) ipoint(6,it) = ipt(2*ix-1,2*iy) endif if (iy .ne. 0 .and. ix .ne. 2**n) then it = it + 1 ipoint(1,it) = ipt(2*ix,2*iy) ipoint(2,it) = ipt(2*ix,2*iy-2) ipoint(3,it) = ipt(2*ix+2,2*iy) ipoint(4,it) = ipt(2*ix,2*iy-1) ipoint(5,it) = ipt(2*ix+1,2*iy-1) ipoint(6,it) = ipt(2*ix+1,2*iy) endif enddo enddo contains function ipt(ix,iy) ! Triangle indexing of (ix,iy) for ix >= iy ! Indexing starts from (0,0) --> 1. implicit none integer*4 :: ipt,ix,iy ipt = (ix*(ix+1))/2 + iy + 1 return end function ipt end ! -------------------------------------------------------------------------- subroutine tri_int(e,q,v,iopt) ! ! Input integration boundary is given by the expansion based on the ! input fermi energy e and the energy profile coefficients q(1..6): ! ! ee(x,y) = q(1) + q(2)*x + q(3)*y + q(4)*x**2 + q(5)*x*y + q(6)*y**2 = e ! ! This curve is a conic section or degenerate reduction thereof ! ! Find affine transformation that places conic in standard position ! ! Input simplex has vertices at A=(0,0), B=(0,1), and C=(1,0). ! Apply affine transformation to simplex and its intersections with ! boundary curve ! ! Identify portions of simplex and boundary curve surrounding the ! occupied portion of the brillouin zone. ! ! Using line-integral methods evaluate integrals (in standard conic ! coordinates) of 1, x, y, x**2, x*y, y**2 ! ! Apply inverse of the affine transformation--transformed integrals ! are output v(1..6) ! ! Cases controlled by iopt: ! iopt= 1 Normal operation; output v(1..6) ! iopt= 0 Density of states; output v(1) only ! iopt=-1 Compute only v(1) (simplex total occupancy) ! ! ---------------------------------------------------------------------- ! implicit none real*8 :: e,q,v,vt,eps,zero,one,two,four,third,qtr,half real*8 :: fifth,sevth,q1,q2,q3,q4,q6,p1,p2,qq,cc,ss,t real*8 :: b1,b2,x1,x2,y1,y2,v1,v2,v4 integer*4 :: iopt,ilst,max_s,iq,i,j,max_q,k,max_i,ix dimension q(6),v(6) dimension vt(6),ilst(3,12) logical isx parameter (eps=1.d-10) ! provisional parameter (zero=0.d0,half=0.5d0,one=1.d0,two=2.d0,four=4.d0) parameter (third=1.d0/3.d0,qtr=1.d0/4.d0,fifth=1.d0/5.d0) parameter (sevth=1.d0/7.d0) type ptlist real*8 :: x real*8 :: y integer*4 :: id ! sequence numbers logical :: tst ! "used" flag in integrations end type ptlist type (ptlist) spt(12),qpt(12),tpt max_s = 1 ! counter for intersection points q1 = q(1) - e if (abs(q1) < eps) q1 = -eps ! Avoid intersections at vertices q2 = q(2); q3 = q(3) if (abs(q1+q3+q(6)) < eps) q3 = -q1-q(6)-eps if (abs(q1+q2+q(4)) < eps) q2 = -q1-q(4)-eps isx = .false. ! Switch for line integration on if (q1 < zero) isx = .true. ! simplex (T if vertex A below e) spt%id = (/ (i, i=1,12) /) spt(1)%x = zero ! Vertex A onto critical-point list spt(1)%y = zero spt(1)%id = 0 ! Find points from A to B (in order) where ee=e: call edgept(q1,q3,q(6)) if (p1 .ne. zero) then max_s = max_s + 1 spt(max_s)%x = zero spt(max_s)%y = p1 if (p2 .ne. zero) then max_s = max_s + 1 spt(max_s)%x = zero spt(max_s)%y = p2 endif endif max_s = max_s + 1 spt(max_s)%x = zero ! Add vertex B to list spt(max_s)%y = one spt(max_s)%id = 0 ! Next along B-C: call edgept(q1+q3+q(6),q2-q3+q(5)-two*q(6),q(4)-q(5)+q(6)) if (p1 .ne. zero) then max_s = max_s + 1 spt(max_s)%x = p1 spt(max_s)%y = one - p1 if (p2 .ne. zero) then max_s = max_s + 1 spt(max_s)%x = p2 spt(max_s)%y = one - p2 endif endif max_s = max_s + 1 ! Add vertex C to list spt(max_s)%x = one spt(max_s)%y = zero spt(max_s)%id = 0 ! Finally along C-A: call edgept(q1,q2,q(4)) if (p2 .ne. 0) then max_s = max_s + 1 spt(max_s)%x = p2 spt(max_s)%y = zero endif if (p1 .ne. 0) then max_s = max_s + 1 spt(max_s)%x = p1 spt(max_s)%y = zero endif ! Bring conic section to standard form and assign identifier if (abs(q(5)) < eps) then if (abs(q(6)) > abs(q(4))) then ! either no rotation or 90 deg. q4 = q(6) q6 = q(4) cc = zero ss = one else q4 = q(4) q6 = q(6) cc = one ss = zero endif else ! rotate to remove x*y term qq = q(4) + q(6) ! (making |q4| >= |q6|) q4 = half*(qq+sign(sqrt((q(4)-q(6))**2+q(5)**2),qq)) q6 = qq - q4 t = two*(q4 - q(4))/q(5) cc = one/sqrt(one + t**2) ss = t*cc endif q2 = cc*q(2) + ss*q(3) ! rotate linear terms q3 = cc*q(3) - ss*q(2) if (abs(q4) < eps) then b1 = zero ! No quadratic terms. Origin b2 = zero ! shift not performed qq = sqrt(q(2)**2 + q(3)**2) if (qq >= eps) then cc = q(2)/qq ss = q(3)/qq ! computed rotation iq = 5 ! vertical line(s) else go to 800 ! flat surface--special coding endif else b1 = half*q2/q4 ! Origin shift: remove x term q1 = q1 - half*q2*b1 if (abs(q6) >= eps) then b2 = half*q3/q6 ! Origin shift: remove y term q1 = q1 - half*q3*b2 if (q1 == zero) q1 = eps*q4 ! Avoid some degeneracies if (q4*q6 > zero) then if (q1*q4 < zero) then iq = 1 ! ellipse else iq = 0 endif else if (q1*q4 >= zero) then iq = 2 ! hyperbola (foci on y axis) else iq = 3 ! hyperbola (foci on x axis) endif else if (abs(q3) >= eps) then b2 = (q1 - q4)/q3 ! Force origin outside of curve q1 = q4 ! (aids in finding path direction) iq = 4 ! parabola (directrix || x axis) else b2 = zero iq = 5 ! vertical line(s) endif endif endif ! Transform simplex critical/intersection points to standard conic coordinates do j = 1,max_s p1 = cc*spt(j)%x + ss*spt(j)%y p2 = cc*spt(j)%y - ss*spt(j)%x spt(j)%x = p1 + b1 spt(j)%y = p2 + b2 enddo ! Make list of conic intersection points max_q = 0 do j = 1,max_s if (spt(j)%id == 0) cycle max_q = max_q + 1 qpt(max_q) = spt(j) enddo ! Add conic critical points (if any) if (iq == 1 .or. iq == 3) then qq = sqrt(-q1/q4) max_q = max_q + 2 qpt(max_q-1)%x = qq qpt(max_q-1)%y = zero qpt(max_q-1)%id = 0 qpt(max_q)%x = -qq qpt(max_q)%y = zero qpt(max_q)%id = 0 endif spt%id = 0 ! Will reload with sequencing id's ! Sort conic critical/intersection points into path order do j = 1,max_q ! Need j=max_q for sequencing line do i = j+1,max_q if (iq == 3 .or. iq == 5) then if (qpt(j)%x >= zero .and. & & (qpt(i)%x <= zero .or. qpt(i)%y > qpt(j)%y) & & .or. qpt(i)%x <= zero .and. qpt(i)%y < qpt(j)%y) cycle else if (qpt(j)%y >= zero .and. & & (qpt(i)%y < zero .or. qpt(i)%x < qpt(j)%x) & & .or. qpt(i)%y < zero .and. qpt(i)%x > qpt(j)%x) cycle endif tpt = qpt(j) qpt(j) = qpt(i) qpt(i) = tpt enddo if (qpt(j)%id /= 0) spt(qpt(j)%id)%id = j ! Set up seq no. in qpt enddo ! Identify integrals over path segments in array ilst k = 1 ! qpt path direction (- if if (q1 < zero) k = -1 ! conic origin below e) spt%tst = .false. ! Mark as unused max_i = 0 ! integral counter i = 1 101 if (spt(i)%tst) go to 105 ! spt point already used--exit 102 spt(i)%tst = .true. ! mark as "used" j = i i = mod(i,max_s) + 1 ! next i (range 1 thru max_s) if (isx) call add_integ(j,i,0) ! edge segment to path list if (spt(i)%id == 0) go to 101 ! spt loop if not root point if (.not. isx) then isx = .true. ! spt loop if root is start go to 101 ! (not end) of side boundary endif i = spt(i)%id ! get qpt corr. to i 103 j = i i = mod(i+k+max_q-1,max_q) + 1 ! next i ix = 1 ! +1, upper branch; -1, lower if (qpt(i)%y < zero .or. qpt(j)%y < 0 .or. qpt(i)%y==0 & & .and. qpt(j)%y==0 .and. k*fp(j,i,1) > zero) ix = -1 call add_integ(j,i,ix) ! curve segment to path list if (qpt(i)%id == 0) go to 103 ! qpt loop if not root point i = qpt(i)%id ! get spt corr. to i go to 101 ! next spt point 105 if (max_s == 3) then if (iq==1) then x1 = qq - b1 y1 = -b2 x2 = cc*x1 - ss*y1 y2 = cc*y1 + ss*x1 if (x2 > zero .and. x2 < one .and. y2 > zero .and. y2 < one & & .and. x2+y2 < one) then ix = -1; if (isx) ix=1 call add_integ(1,2,ix) ! integrate over interior call add_integ(2,1,-ix) ! ellipse endif endif else isx = .false. do i = 1,max_s if (.not. spt(i)%tst) go to 102 ! Find and process any enddo ! unused spt points endif ! Evaluate integrals indicated by array ilst; if iopt<0 only area--vt(1) ! if iopt=0 only DOS vt = zero ! initialize if (iq <= 3) then q1 = -q1/q6 ! modified defs for conic q4 = -q4/q6 ! integrals else if (iq == 4) then q1 = -q1/q3 q4 = -q4/q3 endif do i = 1,max_i j = ilst(1,i) k = ilst(2,i) ix = ilst(3,i) ! 0 for straight line, +/-1 for conics if (iopt == 0) then ! DOS case if (ix == 0) cycle if (iq <= 3) then p1 = ix*half/(q6*sqrt(abs(q4))) if (iq == 1) then v1 = p1*(asin(f1(k))-asin(f1(j))) else v1 = p1*log((f1(k)+f2(k))/(f1(j)+f2(j))) endif vt(1) = vt(1) + v1 v2 = half*ix*(fy(k)-fy(j))/q6 vt(2) = vt(2) + v2 vt(3) = vt(3) + half*fp(j,k,1)/q6 vt(4) = vt(4) - half*((q1/q4)*v1 - v2) vt(5) = vt(5) + qtr*fp(j,k,2)/q6 vt(6) = vt(6) + half*(q1*v1 + q4*v2) else if (iq == 4) then vt(1) = vt(1) + fp(j,k,1)/q3 vt(2) = vt(2) + half*fp(j,k,2)/q3 vt(3) = vt(3) + (third*q4*fp(j,k,3) + q1*fp(j,k,1))/q3 vt(4) = vt(4) + third*fp(j,k,3)/q3 vt(5) = vt(5) + (qtr*q4*fp(j,k,4) + half*q1*fp(j,k,2))/q3 vt(6) = vt(6) + (fifth*q4**2*fp(j,k,5) + & & q1*(two*third*q4*fp(j,k,3) + q1*fp(j,k,1)))/q3 endif else ! Normal case if (ix == 0) then ! straight-line segments if (abs(fs(j,k,1)) .lt. eps) cycle p2 = (spt(k)%y-spt(j)%y)/fs(j,k,1) p1 = spt(j)%y - p2*spt(j)%x v1 = p1*fs(j,k,1) + half*p2*fs(j,k,2) vt(1) = vt(1) + v1 if (iopt > 0) then v2 = half*p1*fs(j,k,2) + third*p2*fs(j,k,3) vt(2) = vt(2) + v2 vt(3) = vt(3) + half*(p2*v2 + p1*v1) v4 = third*p1*fs(j,k,3) + qtr*p2*fs(j,k,4) vt(4) = vt(4) + v4 vt(5) = vt(5) + half*(p2*v4 + p1*v2) vt(6) = vt(6) + third*(p2*(p2*v4+two*p1*v2) + p1**2*v1) endif else if (iq <= 3) then ! ellipse or hyperbola p1 = ix*half*abs(q1)/sqrt(abs(q4)) v1 = f1(k)*f2(k)-f1(j)*f2(j) if (iq == 1 .and. abs(p1) .gt. eps) then v1 = p1*(v1 + asin(f1(k))-asin(f1(j))) else if (abs(p1) .gt. eps) then v1=p1*(v1+log((f1(k)+f2(k))/(f1(j)+f2(j)))*sign(one,q1)) else v1 = 0.d0 endif vt(1) = vt(1) + v1 if (iopt > 0) then vt(2) = vt(2) + ix*third*(fx(k) - fx(j)) vt(3) = vt(3) + half*(third*q4*fp(j,k,3)+q1*fp(j,k,1)) v4 = qtr*(ix*(qpt(k)%x*fx(k)-qpt(j)%x*fx(j)) - q1*v1/q4) vt(4) = vt(4) + v4 vt(5) = vt(5) + qtr*(half*q4*fp(j,k,4) + q1*fp(j,k,2)) vt(6) = vt(6) + third*(q1*v1 + q4*v4) endif else ! parabola, or vert line-->0 v1=third*q4*fp(j,k,3) + q1*fp(j,k,1) vt(1) = vt(1) + v1 if (iopt > 0) then vt(2) = vt(2) + qtr*q4*fp(j,k,4) + half*q1*fp(j,k,2) v4 = fifth*q4*fp(j,k,5) + third*q1*fp(j,k,3) vt(3) = vt(3) + half*(q4*v4 + q1*v1) vt(4) = vt(4) + v4 vt(5) = vt(5) + qtr*(third*q4**2*fp(j,k,6) & & + q1*q4*fp(j,k,4) + q1**2*fp(j,k,2)) vt(6) = vt(6) + q1*q4*v4 + third*(sevth*q4**3*fp(j,k,7) & & + q1**3*fp(j,k,1)) endif endif endif enddo v(1) = vt(1) if (iopt == -1) return ! Transform vt(2)..vt(6) to original coordinate system vt(2) = vt(2) - b1*v(1) vt(3) = vt(3) - b2*v(1) v(2) = cc*vt(2) - ss*vt(3) v(3) = cc*vt(3) + ss*vt(2) vt(4) = vt(4) - two*b1*vt(2) - b1**2*v(1) vt(5) = vt(5) - b1*vt(3) - b2*vt(2) - b1*b2*v(1) vt(6) = vt(6) - two*b2*vt(3) - b2**2*v(1) v(4) = cc**2*vt(4) + ss**2*vt(6) - two*cc*ss*vt(5) v(5) = (cc**2 - ss**2)*vt(5) + cc*ss*(vt(4) - vt(6)) v(6) = vt(4) + vt(6) - v(4) return ! Flat surface (special case) 800 if (q(1) == e .and. iopt == 0) then v = 1.d10 else if (q(1) <= e .and. iopt .ne. 0) then v(1) = half ! Assume full occupancy if (iopt == 1) then v(2) = half*third v(3) = v(2) v(4) = half*v(2) v(5) = half*v(4) v(6) = v(4) endif else v = zero endif contains subroutine edgept(c,b,a) ! Finds roots p1,p2 of a*z**2 + b*z + c such that 0 < root < 1. Ignores ! double roots completely. If no such roots, sets p1=0. If one ! such root, sets p1=root and p2=0. If two roots, they are p1 < p2. implicit none real*8 :: a,b,c,s,t p1 = zero p2 = zero if (a == zero) then if (b == zero) return p1 = -c/b else t = b**2 - four*a*c if (t <= zero) return t = sqrt(t)/abs(a) s = -b/a p1 = half*(s-t) p2 = half*(s+t) if (p2 <= zero .or. p2 >= one) p2 = zero endif if (p1 <= zero .or. p1 >= one) then p1 = p2 p2 = zero endif end subroutine edgept subroutine add_integ(j,i,m) implicit none integer*4 :: j,i,m max_i = max_i + 1 ilst(1,max_i) = j ilst(2,max_i) = i ilst(3,max_i) = m end subroutine add_integ function fx(i) implicit none real*8 :: fx integer*4 :: i fx = max(q1+q4*qpt(i)%x**2,zero)**(1.5)/q4 end function fx function fy(i) implicit none real*8 :: fy integer*4 :: i fy = max(q1+q4*qpt(i)%x**2,zero)**(0.5)/q4 end function fy function fp(i,j,k) implicit none real*8 :: fp integer*4 :: i,j,k fp = qpt(j)%x**k-qpt(i)%x**k end function fp function fs(i,j,k) implicit none real*8 :: fs integer*4 :: i,j,k fs = spt(j)%x**k-spt(i)%x**k end function fs function f1(i) implicit none real*8 :: f1 integer*4 :: i f1 = sqrt(abs(q4/q1))*qpt(i)%x if (iq == 1) f1 = sign(min(abs(f1),one),f1) end function f1 function f2(i) implicit none real*8 :: f2 integer*4 :: i f2 = sqrt(max((q4*qpt(i)%x**2+q1)/abs(q1),zero)) end function f2 end