--- /dev/null
+ subroutine integral(gamma1,gamma2,gamma3,gamma4,ity1,ity2,a1,a2,
+ & si1,si2,si3,si4,transp,q)
+ implicit none
+ integer ity1,ity2
+ integer ilam1,ilam2,ilam3,ilam4,iincr
+ double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
+ & lambda2,lambda3,lambda4
+ logical transp
+ double precision elocal,ele
+ double precision delta,delta2,sum,ene,sumene,boltz
+ double precision q,a1(2,2),a2(2,2),si1,si2,si3,si4
+ double precision conv /.01745329252d0/,pi /3.141592654d0/
+
+ iincr=20
+ delta=iincr*conv
+ delta2=0.5d0*delta
+cd print *,'iincr',iincr,' delta',delta
+cd write(2,*) gamma1,gamma2,ity1,ity2,a1,a2,si1,si2,si3,si4,transp
+
+cd do ilam1=-180,180,5
+cd do ilam2=-180,180,5
+cd lambda1=ilam1*conv+delta2
+cd lambda2=ilam2*conv+delta2
+cd write(2,'(2i5,2f10.5)') ilam1,ilam2,elocal(2,lambda1,lambda2),
+cd & ele(lambda1,lambda2,a1,1.0d0,1.d00)
+cd enddo
+cd enddo
+cd stop
+
+ sum=0.0d0
+ sumene=0.0d0
+ do ilam1=-180,179,iincr
+ do ilam2=-180,179,iincr
+ do ilam3=-180,179,iincr
+ do ilam4=-180,179,iincr
+ lambda1=ilam1*conv+delta2
+ lambda2=ilam2*conv+delta2
+ lambda3=ilam3*conv+delta2
+ lambda4=ilam4*conv+delta2
+cd write (2,*) ilam1,ilam2,ilam3,ilam4
+cd write (2,*) lambda1,lambda2,lambda3,lambda4
+ ene=
+ & -elocal(ity1,lambda1,lambda2,.false.)*
+ & elocal(ity2,lambda3,lambda4,transp)*
+ & ele(si1*lambda1+gamma1,si3*lambda3+gamma3,a1)*
+ & ele(si2*lambda2+gamma2,si4*lambda4+gamma4,a2)
+cd write (2,*) elocal(ity1,lambda1,gamma1-pi-lambda2),
+cd & elocal(ity2,lambda3,gamma2-pi-lambda4),
+cd & ele(lambda1,lambda2,a1,si1,si3),
+cd & ele(lambda3,lambda4,a2,si2,si4)
+ sum=sum+ene
+ enddo
+ enddo
+ enddo
+ enddo
+ q=sum/(2*pi)**4*delta**4
+ write (2,* )'sum',sum,' q',q
+ return
+ end
+c---------------------------------------------------------------------------
+ subroutine integral3(gamma1,gamma2,ity1,ity2,ity3,ity4,
+ & a1,koniec,q1,q2,q3,q4)
+ implicit none
+ integer ity1,ity2,ity3,ity4
+ integer ilam1,ilam2,ilam3,ilam4,iincr
+ double precision gamma1,gamma2,gamma3,gamma4,beta,lambda1,
+ & lambda2,lambda3,lambda4
+ logical koniec
+ double precision elocal,ele
+ double precision delta,delta2,sum1,sum2,sum3,sum4,
+ & ene1,ene2,ene3,ene4,boltz
+ double precision q1,q2,q3,q4,a1(2,2),a2(2,2)
+ double precision conv /.01745329252d0/,pi /3.141592654d0/
+
+ iincr=60
+ delta=iincr*conv
+ delta2=0.5d0*delta
+cd print *,'iincr',iincr,' delta',delta
+ write(2,*) gamma1,gamma2,ity1,ity2,ity3,ity4,a1,koniec
+
+cd do ilam1=-180,180,5
+cd do ilam2=-180,180,5
+cd lambda1=ilam1*conv+delta2
+cd lambda2=ilam2*conv+delta2
+cd write(2,'(2i5,2f10.5)') ilam1,ilam2,elocal(2,lambda1,lambda2),
+cd & ele(lambda1,lambda2,a1,1.0d0,1.d00)
+cd enddo
+cd enddo
+cd stop
+
+ sum1=0.0d0
+ sum2=0.0d0
+ sum3=0.0d0
+ sum4=0.0d0
+ do ilam1=-180,179,iincr
+ do ilam2=-180,179,iincr
+ do ilam3=-180,179,iincr
+ do ilam4=-180,179,iincr
+ lambda1=ilam1*conv+delta2
+ lambda2=ilam2*conv+delta2
+ lambda3=ilam3*conv+delta2
+ lambda4=ilam4*conv+delta2
+cd write (2,*) ilam1,ilam2,ilam3,ilam4
+cd write (2,*) lambda1,lambda2,lambda3,lambda4
+ if (.not.koniec) then
+ ene1=
+ & elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
+ & elocal(ity3,lambda3,gamma2-pi-lambda4,.false.)*
+ & ele(lambda2,lambda4,a1)
+ else
+ ene1=
+ & elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
+ & elocal(ity3,lambda3,lambda4,.false.)*
+ & ele(lambda2,-lambda4,a1)
+ endif
+ ene2=
+ & elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
+ & elocal(ity4,lambda3,lambda4,.false.)*
+ & ele(lambda2,lambda3,a1)
+ if (.not.koniec) then
+ ene3=
+ & elocal(ity2,lambda1,lambda2,.false.)*
+ & elocal(ity3,lambda3,gamma2-pi-lambda4,.false.)*
+ & ele(lambda1,lambda4,a1)
+ else
+ ene3=
+ & elocal(ity2,lambda1,lambda2,.false.)*
+ & elocal(ity3,lambda3,lambda4,.false.)*
+ & ele(lambda1,-lambda4,a1)
+ endif
+ ene4=
+ & elocal(ity2,lambda1,lambda2,.false.)*
+ & elocal(ity4,lambda3,lambda4,.false.)*
+ & ele(lambda1,lambda3,a1)
+ sum1=sum1+ene1
+ sum2=sum2+ene2
+ sum3=sum3+ene3
+ sum4=sum4+ene4
+ enddo
+ enddo
+ enddo
+ enddo
+ q1=sum1/(2*pi)**4*delta**4
+ q2=sum2/(2*pi)**4*delta**4
+ q3=sum3/(2*pi)**4*delta**4
+ q4=sum4/(2*pi)**4*delta**4
+ write (2,* )'sum',sum1,sum2,sum3,sum4,' q',q1,q2,q3,q4
+ return
+ end
+c-------------------------------------------------------------------------
+ subroutine integral5(gamma1,gamma2,gamma3,gamma4,ity1,ity2,ity3,
+ & ity4,ity5,ity6,a1,a2,si1,si2,si3,si4,transp,ene1,ene2,ene3,ene4)
+ implicit none
+ integer ity1,ity2,ity3,ity4,ity5,ity6
+ integer ilam1,ilam2,ilam3,ilam4,ilam5,iincr
+ double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
+ & lambda2,lambda3,lambda4,lambda5
+ logical transp
+ double precision elocal,ele
+ double precision eloc1,eloc2,eloc3,eloc4,eloc5,eloc6,ele1,ele2
+ double precision delta,delta2,sum,ene,sumene,pom
+ double precision ene1,ene2,ene3,ene4,sum1,sum2,sum3,sum4,
+ & a1(2,2),a2(2,2)
+ integer si1,si2,si3,si4
+ double precision conv /.01745329252d0/,pi /3.141592654d0/
+
+ iincr=60
+ delta=iincr*conv
+ delta2=0.5d0*delta
+cd print *,'iincr',iincr,' delta',delta
+cd write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
+cd & ' gamma3=',gamma3,' gamma4=',gamma4
+cd write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
+cd write(2,*) 'a1=',a1
+cd write(2,*) 'a2=',a2
+cd write(2,*) si1,si2,si3,si4,transp
+
+ sum1=0.0d0
+ sum2=0.0d0
+ sum3=0.0d0
+ sum4=0.0d0
+ do ilam1=-180,179,iincr
+ do ilam2=-180,179,iincr
+ do ilam3=-180,179,iincr
+ do ilam4=-180,179,iincr
+ do ilam5=-180,179,iincr
+ lambda1=ilam1*conv+delta2
+ lambda2=ilam2*conv+delta2
+ lambda3=ilam3*conv+delta2
+ lambda4=ilam4*conv+delta2
+ lambda5=ilam5*conv+delta2
+ if (transp) then
+ ele1=ele(lambda1,si4*lambda4,a1)
+ ele2=ele(lambda2,lambda3,a2)
+ else
+ ele1=ele(lambda1,lambda3,a1)
+ ele2=ele(lambda2,si4*lambda4,a2)
+ endif
+ eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
+ eloc5=elocal(ity5,lambda3,gamma4-pi-si4*lambda4,.false.)
+ pom=ele1*ele2*eloc2*eloc5
+ if (si1.gt.0) then
+ eloc1=elocal(ity1,lambda5,gamma1-pi-lambda1,.false.)
+ sum1=sum1+pom*eloc1
+ endif
+ eloc3=elocal(ity3,lambda2,lambda5,.false.)
+ sum2=sum2+pom*eloc3
+ eloc4=elocal(ity4,lambda5,gamma3-pi-lambda3,.false.)
+ sum3=sum3+pom*eloc4
+ if (si4.gt.0) then
+ eloc6=elocal(ity6,lambda4,lambda5,.false.)
+ sum4=sum4+pom*eloc6
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ pom=1.0d0/(2*pi)**5*delta**5
+ ene1=sum1*pom
+ ene2=sum2*pom
+ ene3=sum3*pom
+ ene4=sum4*pom
+c write (2,* )'sum',sum1,sum2,sum3,sum4,' q',ene1,ene2,ene3,ene4
+ return
+ end
+c-------------------------------------------------------------------------
+ subroutine integral_turn6(gamma1,gamma2,gamma3,gamma4,ity1,ity2,
+ & ity3,ity4,ity5,ity6,a1,a2,ene_turn6)
+ implicit none
+ integer ity1,ity2,ity3,ity4,ity5,ity6
+ integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
+ double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
+ & lambda2,lambda3,lambda4,lambda5,lambda6
+ logical transp
+ double precision elocal,ele
+ double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
+ & eloc61,ele1,ele2
+ double precision delta,delta2,sum,ene,sumene,pom,ene5
+ double precision ene_turn6,sum5,a1(2,2),a2(2,2)
+ double precision conv /.01745329252d0/,pi /3.141592654d0/
+
+ iincr=60
+ delta=iincr*conv
+ delta2=0.5d0*delta
+cd print *,'iincr',iincr,' delta',delta
+ write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
+ & ' gamma3=',gamma3,' gamma4=',gamma4
+ write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
+ write(2,*) 'a1=',a1
+ write(2,*) 'a2=',a2
+
+ sum5=0.0d0
+ do ilam1=-180,179,iincr
+ do ilam2=-180,179,iincr
+ do ilam3=-180,179,iincr
+ do ilam4=-180,179,iincr
+ do ilam5=-180,179,iincr
+ lambda1=ilam1*conv+delta2
+ lambda2=ilam2*conv+delta2
+ lambda3=ilam3*conv+delta2
+ lambda4=ilam4*conv+delta2
+ lambda5=ilam5*conv+delta2
+ ele1=ele(lambda1,-lambda4,a1)
+ ele2=ele(lambda2,lambda3,a2)
+ eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
+ eloc5=elocal(ity5,lambda3,lambda4,.false.)
+ pom=ele1*ele2*eloc2*eloc5
+ eloc3=elocal(ity3,lambda2,gamma3-pi-lambda5,.false.)
+ eloc4=elocal(ity4,lambda5,gamma4-pi-lambda3,.false.)
+ sum5=sum5+pom*eloc3*eloc4
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ pom=-1.0d0/(2*pi)**5*delta**5
+ ene_turn6=sum5*pom
+c print *,'sum6',sum6,' ene6',ene6
+ return
+ end
+c-------------------------------------------------------------------------
+ subroutine integral6(gamma1,gamma2,gamma3,gamma4,ity1,ity2,ity3,
+ & ity4,ity5,ity6,a1,a2,si1,si2,si3,si4,transp,ene1,ene2,ene3,ene4,
+ & ene5,ene6)
+ implicit none
+ integer ity1,ity2,ity3,ity4,ity5,ity6
+ integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
+ double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
+ & lambda2,lambda3,lambda4,lambda5,lambda6
+ logical transp
+ double precision elocal,ele
+ double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
+ & eloc61,ele1,ele2
+ double precision delta,delta2,sum,ene,sumene,pom
+ double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
+ & sum4,sum5,sum6,a1(2,2),a2(2,2)
+ integer si1,si2,si3,si4
+ double precision conv /.01745329252d0/,pi /3.141592654d0/
+
+ iincr=60
+ delta=iincr*conv
+ delta2=0.5d0*delta
+cd print *,'iincr',iincr,' delta',delta
+cd write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
+cd & ' gamma3=',gamma3,' gamma4=',gamma4
+cd write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
+cd write(2,*) 'a1=',a1
+cd write(2,*) 'a2=',a2
+cd write(2,*) si1,si2,si3,si4,transp
+
+ sum1=0.0d0
+ sum2=0.0d0
+ sum3=0.0d0
+ sum4=0.0d0
+ sum5=0.0d0
+ sum6=0.0d0
+ eloc1=0.0d0
+ eloc6=0.0d0
+ eloc61=0.0d0
+ do ilam1=-180,179,iincr
+ do ilam2=-180,179,iincr
+ do ilam3=-180,179,iincr
+ do ilam4=-180,179,iincr
+ do ilam5=-180,179,iincr
+ do ilam6=-180,179,iincr
+ lambda1=ilam1*conv+delta2
+ lambda2=ilam2*conv+delta2
+ lambda3=ilam3*conv+delta2
+ lambda4=ilam4*conv+delta2
+ lambda5=ilam5*conv+delta2
+ lambda6=ilam6*conv+delta2
+ if (transp) then
+ ele1=ele(lambda1,si4*lambda4,a1)
+ ele2=ele(lambda2,lambda3,a2)
+ else
+ ele1=ele(lambda1,lambda3,a1)
+ ele2=ele(lambda2,si4*lambda4,a2)
+ endif
+ eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
+ eloc5=elocal(ity5,lambda3,gamma4-pi-si4*lambda4,.false.)
+ pom=ele1*ele2*eloc2*eloc5
+ if (si1.gt.0) then
+ eloc1=elocal(ity1,lambda5,gamma1-pi-lambda1,.false.)
+ endif
+ eloc3=elocal(ity3,lambda2,lambda6,.false.)
+ sum1=sum1+pom*eloc1*eloc3
+ eloc4=elocal(ity4,lambda5,gamma3-pi-lambda3,.false.)
+ if (si4.gt.0) then
+ eloc6=elocal(ity6,lambda4,lambda6,.false.)
+ eloc61=elocal(ity6,lambda4,lambda5,.false.)
+ endif
+ sum2=sum2+pom*eloc4*eloc6
+ eloc41=elocal(ity4,lambda6,gamma3-pi-lambda3,.false.)
+ sum3=sum3+pom*eloc1*eloc41
+ sum4=sum4+pom*eloc1*eloc6
+ sum5=sum5+pom*eloc3*eloc4
+ sum6=sum6+pom*eloc3*eloc61
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ pom=-1.0d0/(2*pi)**6*delta**6
+ ene1=sum1*pom
+ ene2=sum2*pom
+ ene3=sum3*pom
+ ene4=sum4*pom
+ ene5=sum5*pom
+ ene6=sum6*pom
+c print *,'sum6',sum6,' ene6',ene6
+ return
+ end
+c-------------------------------------------------------------------------
+ subroutine integral3a(gamma1,gamma2,ity1,ity2,a1,si1,ene1)
+ implicit none
+ integer ity1,ity2,ity3,ity4,ity5,ity6
+ integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
+ double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
+ & lambda2,lambda3,lambda4,lambda5,lambda6
+ logical transp
+ double precision elocal,ele
+ double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
+ & eloc61,ele1,ele2
+ double precision delta,delta2,sum,ene,sumene,pom
+ double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
+ & sum4,sum5,sum6,a1(2,2),a2(2,2)
+ integer si1,si2,si3,si4
+ double precision conv /.01745329252d0/,pi /3.141592654d0/
+
+ iincr=60
+ delta=iincr*conv
+ delta2=0.5d0*delta
+cd print *,'iincr',iincr,' delta',delta
+cd write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2
+cd write(2,*) ity1,ity2
+cd write(2,*) 'a1=',a1
+cd write(2,*) si1,
+
+ sum1=0.0d0
+ eloc1=0.0d0
+ do ilam1=-180,179,iincr
+ do ilam2=-180,179,iincr
+ do ilam3=-180,179,iincr
+ lambda1=ilam1*conv+delta2
+ lambda2=ilam2*conv+delta2
+ lambda3=ilam3*conv+delta2
+ ele1=ele(lambda1,si1*lambda3,a1)
+ eloc1=elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)
+ if (si1.gt.0) then
+ eloc2=elocal(ity2,lambda2,gamma2-pi-lambda3,.false.)
+ else
+ eloc2=elocal(ity2,lambda2,lambda3,.false.)
+ endif
+ sum1=sum1+ele1*eloc1*eloc2
+ enddo
+ enddo
+ enddo
+ pom=1.0d0/(2*pi)**3*delta**3
+ ene1=sum1*pom
+ return
+ end
+c-------------------------------------------------------------------------
+ subroutine integral4a(gamma1,gamma2,gamma3,ity1,ity2,ity3,a1,si1,
+ & ene1)
+ implicit none
+ integer ity1,ity2,ity3,ity4,ity5,ity6
+ integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
+ double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
+ & lambda2,lambda3,lambda4,lambda5,lambda6
+ logical transp
+ double precision elocal,ele
+ double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
+ & eloc61,ele1,ele2
+ double precision delta,delta2,sum,ene,sumene,pom
+ double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
+ & sum4,sum5,sum6,a1(2,2),a2(2,2)
+ integer si1,si2,si3,si4
+ double precision conv /.01745329252d0/,pi /3.141592654d0/
+
+ iincr=60
+ delta=iincr*conv
+ delta2=0.5d0*delta
+cd print *,'iincr',iincr,' delta',delta
+cd write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
+cd & ' gamma3=',gamma3
+cd write(2,*) ity1,ity2,ity3
+cd write(2,*) 'a1=',a1
+cd write(2,*) 'si1=',si1
+ sum1=0.0d0
+ do ilam1=-180,179,iincr
+ do ilam2=-180,179,iincr
+ do ilam3=-180,179,iincr
+ do ilam4=-180,179,iincr
+ lambda1=ilam1*conv+delta2
+ lambda2=ilam2*conv+delta2
+ lambda3=ilam3*conv+delta2
+ lambda4=ilam4*conv+delta2
+ ele1=ele(lambda1,si1*lambda4,a1)
+ eloc1=elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)
+ eloc2=elocal(ity2,lambda2,gamma2-pi-lambda3,.false.)
+ if (si1.gt.0) then
+ eloc3=elocal(ity3,lambda3,gamma3-pi-lambda4,.false.)
+ else
+ eloc3=elocal(ity3,lambda3,lambda4,.false.)
+ endif
+ sum1=sum1+ele1*eloc1*eloc2*eloc3
+ enddo
+ enddo
+ enddo
+ enddo
+ pom=-1.0d0/(2*pi)**4*delta**4
+ ene1=sum1*pom
+ return
+ end
+c-------------------------------------------------------------------------
+ double precision function elocal(i,x,y,transp)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.TORSION'
+ integer i
+ double precision x,y,u(2),v(2),cu(2),dv(2),ev(2)
+ double precision scalar2
+ logical transp
+ u(1)=dcos(x)
+ u(2)=dsin(x)
+ v(1)=dcos(y)
+ v(2)=dsin(y)
+ if (transp) then
+ call matvec2(cc(1,1,i),v,cu)
+ call matvec2(dd(1,1,i),u,dv)
+ call matvec2(ee(1,1,i),u,ev)
+ elocal=scalar2(b1(1,i),v)+scalar2(b2(1,i),u)+scalar2(cu,v)+
+ & scalar2(dv,u)+scalar2(ev,v)
+ else
+ call matvec2(cc(1,1,i),u,cu)
+ call matvec2(dd(1,1,i),v,dv)
+ call matvec2(ee(1,1,i),v,ev)
+ elocal=scalar2(b1(1,i),u)+scalar2(b2(1,i),v)+scalar2(cu,u)+
+ & scalar2(dv,v)+scalar2(ev,u)
+ endif
+ return
+ end
+c-------------------------------------------------------------------------
+ double precision function ele(x,y,a)
+ implicit none
+ double precision x,y,a(2,2),si1,si2,u(2),v(2),av(2)
+ double precision scalar2
+ u(1)=-cos(x)
+ u(2)= sin(x)
+ v(1)=-cos(y)
+ v(2)= sin(y)
+ call matvec2(a,v,av)
+ ele=scalar2(u,av)
+ return
+ end
projects / unres.git / blobdiff