X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-NEWSC%2Fintlocal.f;fp=source%2Funres%2Fsrc_MD-NEWSC%2Fintlocal.f;h=2dbcc88c1293f9e71736f3f36d730e26f410a035;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_MD-NEWSC/intlocal.f b/source/unres/src_MD-NEWSC/intlocal.f new file mode 100644 index 0000000..2dbcc88 --- /dev/null +++ b/source/unres/src_MD-NEWSC/intlocal.f @@ -0,0 +1,517 @@ + 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