X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_CSA_DiL%2Ffitsq.f;fp=source%2Funres%2Fsrc_CSA_DiL%2Ffitsq.f;h=0000000000000000000000000000000000000000;hp=36cbd3072b4ba678cc567766f1083e417f73b7a0;hb=2a226bfc86eabc6e4eae0c3ad1cbc3cb5417a05a;hpb=a0e685f844163003749ba91dfbf4644bcc8cfa30 diff --git a/source/unres/src_CSA_DiL/fitsq.f b/source/unres/src_CSA_DiL/fitsq.f deleted file mode 100644 index 36cbd30..0000000 --- a/source/unres/src_CSA_DiL/fitsq.f +++ /dev/null @@ -1,364 +0,0 @@ - subroutine fitsq(rms,x,y,nn,t,b,non_conv) - implicit real*8 (a-h,o-z) - include 'COMMON.IOUNITS' -c x and y are the vectors of coordinates (dimensioned (3,n)) of the two -c structures to be superimposed. nn is 3*n, where n is the number of -c points. t and b are respectively the translation vector and the -c rotation matrix that transforms the second set of coordinates to the -c frame of the first set. -c eta = machine-specific variable - - dimension x(3*nn),y(3*nn),t(3) - dimension b(3,3),q(3,3),r(3,3),v(3),xav(3),yav(3),e(3),c(3,3) - logical non_conv -c eta = z00100000 -c small=25.0*rmdcon(3) -c small=25.0*eta -c small=25.0*10.e-10 -c the following is a very lenient value for 'small' - small = 0.0001D0 - non_conv=.false. - fn=nn - do 10 i=1,3 - xav(i)=0.0D0 - yav(i)=0.0D0 - do 10 j=1,3 - 10 b(j,i)=0.0D0 - nc=0 -c - do 30 n=1,nn - do 20 i=1,3 -c write(iout,*)'x = ',x(nc+i),' y = ',y(nc+i) - xav(i)=xav(i)+x(nc+i)/fn - 20 yav(i)=yav(i)+y(nc+i)/fn - 30 nc=nc+3 -c - do i=1,3 - t(i)=yav(i)-xav(i) - enddo - - rms=0.0d0 - do n=1,nn - do i=1,3 - rms=rms+(y(3*(n-1)+i)-x(3*(n-1)+i)-t(i))**2 - enddo - enddo - rms=dabs(rms/fn) - -c write(iout,*)'xav = ',(xav(j),j=1,3) -c write(iout,*)'yav = ',(yav(j),j=1,3) -c write(iout,*)'t = ',(t(j),j=1,3) -c write(iout,*)'rms=',rms - if (rms.lt.small) return - - - nc=0 - rms=0.0D0 - do 50 n=1,nn - do 40 i=1,3 - rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn - do 40 j=1,3 - b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn - 40 c(j,i)=b(j,i) - 50 nc=nc+3 - call sivade(b,q,r,d,non_conv) - sn3=dsign(1.0d0,d) - do 120 i=1,3 - do 120 j=1,3 - 120 b(j,i)=-q(j,1)*r(i,1)-q(j,2)*r(i,2)-sn3*q(j,3)*r(i,3) - call mvvad(b,xav,yav,t) - do 130 i=1,3 - do 130 j=1,3 - rms=rms+2.0*c(j,i)*b(j,i) - 130 b(j,i)=-b(j,i) - if (dabs(rms).gt.small) go to 140 -* write (6,301) - return - 140 if (rms.gt.0.0d0) go to 150 -c write (iout,303) rms - rms=0.0d0 -* stop -c 150 write (iout,302) dsqrt(rms) - 150 continue - return - 301 format (5x,'rms deviation negligible') - 302 format (5x,'rms deviation ',f14.6) - 303 format (//,5x,'negative ms deviation - ',f14.6) - end -c - subroutine sivade(x,q,r,dt,non_conv) - implicit real*8(a-h,o-z) -c computes q,e and r such that q(t)xr = diag(e) - dimension x(3,3),q(3,3),r(3,3),e(3) - dimension h(3,3),p(3,3),u(3,3),d(3) - logical non_conv -c eta = z00100000 -c write (2,*) "SIVADE" - nit = 0 - small=25.0*10.d-10 -c small=25.0*eta -c small=2.0*rmdcon(3) - xnrm=0.0d0 - do 20 i=1,3 - do 10 j=1,3 - xnrm=xnrm+x(j,i)*x(j,i) - u(j,i)=0.0d0 - r(j,i)=0.0d0 - 10 h(j,i)=0.0d0 - u(i,i)=1.0 - 20 r(i,i)=1.0 - xnrm=dsqrt(xnrm) - do 110 n=1,2 - xmax=0.0d0 - do 30 j=n,3 - 30 if (dabs(x(j,n)).gt.xmax) xmax=dabs(x(j,n)) - a=0.0d0 - do 40 j=n,3 - h(j,n)=x(j,n)/xmax - 40 a=a+h(j,n)*h(j,n) - a=dsqrt(a) - den=a*(a+dabs(h(n,n))) - d(n)=1.0/den - h(n,n)=h(n,n)+dsign(a,h(n,n)) - do 70 i=n,3 - s=0.0d0 - do 50 j=n,3 - 50 s=s+h(j,n)*x(j,i) - s=d(n)*s - do 60 j=n,3 - 60 x(j,i)=x(j,i)-s*h(j,n) - 70 continue - if (n.gt.1) go to 110 - xmax=dmax1(dabs(x(1,2)),dabs(x(1,3))) - h(2,3)=x(1,2)/xmax - h(3,3)=x(1,3)/xmax - a=dsqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3)) - den=a*(a+dabs(h(2,3))) - d(3)=1.0/den - h(2,3)=h(2,3)+sign(a,h(2,3)) - do 100 i=1,3 - s=0.0d0 - do 80 j=2,3 - 80 s=s+h(j,3)*x(i,j) - s=d(3)*s - do 90 j=2,3 - 90 x(i,j)=x(i,j)-s*h(j,3) - 100 continue - 110 continue - do 130 i=1,3 - do 120 j=1,3 - 120 p(j,i)=-d(1)*h(j,1)*h(i,1) - 130 p(i,i)=1.0+p(i,i) - do 140 i=2,3 - do 140 j=2,3 - u(j,i)=u(j,i)-d(2)*h(j,2)*h(i,2) - 140 r(j,i)=r(j,i)-d(3)*h(j,3)*h(i,3) - call mmmul(p,u,q) - 150 np=1 - nq=1 - nit=nit+1 -c write (2,*) "nit",nit," e",(x(i,i),i=1,3) - if (nit.gt.10000) then - print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!' - non_conv=.true. - return - endif - if (dabs(x(2,3)).gt.small*(dabs(x(2,2))+abs(x(3,3)))) go to 160 - x(2,3)=0.0d0 - nq=nq+1 - 160 if (dabs(x(1,2)).gt.small*(dabs(x(1,1))+dabs(x(2,2)))) go to 180 - x(1,2)=0.0d0 - if (x(2,3).ne.0.0d0) go to 170 - nq=nq+1 - go to 180 - 170 np=np+1 - 180 if (nq.eq.3) go to 310 - npq=4-np-nq -c write (2,*) "np",np," npq",npq - if (np.gt.npq) go to 230 - n0=0 - do 220 n=np,npq - nn=n+np-1 -c write (2,*) "nn",nn - if (dabs(x(nn,nn)).gt.small*xnrm) go to 220 - x(nn,nn)=0.0d0 - if (x(nn,nn+1).eq.0.0d0) go to 220 - n0=n0+1 -c write (2,*) "nn",nn - go to (190,210,220),nn - 190 do 200 j=2,3 - 200 call givns(x,q,1,j) - go to 220 - 210 call givns(x,q,2,3) - 220 continue -c write (2,*) "nn",nn," np",np," nq",nq," n0",n0 -c write (2,*) "x",(x(i,i),i=1,3) - if (n0.ne.0) go to 150 - 230 nn=3-nq - a=x(nn,nn)*x(nn,nn) - if (nn.gt.1) a=a+x(nn-1,nn)*x(nn-1,nn) - b=x(nn+1,nn+1)*x(nn+1,nn+1)+x(nn,nn+1)*x(nn,nn+1) - c=x(nn,nn)*x(nn,nn+1) - dd=0.5*(a-b) - xn2=c*c - rt=b-xn2/(dd+sign(dsqrt(dd*dd+xn2),dd)) - y=x(np,np)*x(np,np)-rt - z=x(np,np)*x(np,np+1) - do 300 n=np,nn -c write (2,*) "n",n," a",a," b",b," c",c," y",y," z",z - if (dabs(y).lt.dabs(z)) go to 240 - t=z/y - c=1.0/dsqrt(1.0d0+t*t) - s=c*t - go to 250 - 240 t=y/z - s=1.0/dsqrt(1.0d0+t*t) - c=s*t - 250 do 260 j=1,3 - v=x(j,n) - w=x(j,n+1) - x(j,n)=c*v+s*w - x(j,n+1)=-s*v+c*w - a=r(j,n) - b=r(j,n+1) - r(j,n)=c*a+s*b - 260 r(j,n+1)=-s*a+c*b - y=x(n,n) - z=x(n+1,n) - if (dabs(y).lt.dabs(z)) go to 270 - t=z/y - c=1.0/dsqrt(1.0+t*t) - s=c*t - go to 280 - 270 t=y/z - s=1.0/dsqrt(1.0+t*t) - c=s*t - 280 do 290 j=1,3 - v=x(n,j) - w=x(n+1,j) - a=q(j,n) - b=q(j,n+1) - x(n,j)=c*v+s*w - x(n+1,j)=-s*v+c*w - q(j,n)=c*a+s*b - 290 q(j,n+1)=-s*a+c*b - if (n.ge.nn) go to 300 - y=x(n,n+1) - z=x(n,n+2) - 300 continue - go to 150 - 310 do 320 i=1,3 - 320 e(i)=x(i,i) - nit=0 - 330 n0=0 - nit=nit+1 - if (nit.gt.10000) then - print '(a)','!!!! Over 10000 iterations in SIVADE!!!!!' - non_conv=.true. - return - endif -c write (2,*) "e",(e(i),i=1,3) - do 360 i=1,3 - if (e(i).ge.0.0d0) go to 350 - e(i)=-e(i) - do 340 j=1,3 - 340 q(j,i)=-q(j,i) - 350 if (i.eq.1) go to 360 - if (dabs(e(i)).lt.dabs(e(i-1))) go to 360 - call switch(i,1,q,r,e) - n0=n0+1 - 360 continue - if (n0.ne.0) go to 330 -c write (2,*) "e",(e(i),i=1,3) - if (dabs(e(3)).gt.small*xnrm) go to 370 - e(3)=0.0d0 - if (dabs(e(2)).gt.small*xnrm) go to 370 - e(2)=0.0d0 - 370 dt=det(q(1,1),q(1,2),q(1,3))*det(r(1,1),r(1,2),r(1,3)) -c write (2,*) "nit",nit -c write (2,501) (e(i),i=1,3) - return - 501 format (/,5x,'singular values - ',3e15.5) - end - subroutine givns(a,b,m,n) - implicit real*8 (a-h,o-z) - dimension a(3,3),b(3,3) - if (dabs(a(m,n)).lt.dabs(a(n,n))) go to 10 - t=a(n,n)/a(m,n) - s=1.0/dsqrt(1.0+t*t) - c=s*t - go to 20 - 10 t=a(m,n)/a(n,n) - c=1.0/dsqrt(1.0+t*t) - s=c*t - 20 do 30 j=1,3 - v=a(m,j) - w=a(n,j) - x=b(j,m) - y=b(j,n) - a(m,j)=c*v-s*w - a(n,j)=s*v+c*w - b(j,m)=c*x-s*y - 30 b(j,n)=s*x+c*y - return - end - subroutine switch(n,m,u,v,d) - implicit real*8 (a-h,o-z) - dimension u(3,3),v(3,3),d(3) - do 10 i=1,3 - tem=u(i,n) - u(i,n)=u(i,n-1) - u(i,n-1)=tem - if (m.eq.0) go to 10 - tem=v(i,n) - v(i,n)=v(i,n-1) - v(i,n-1)=tem - 10 continue - tem=d(n) - d(n)=d(n-1) - d(n-1)=tem - return - end - subroutine mvvad(b,xav,yav,t) - implicit real*8 (a-h,o-z) - dimension b(3,3),xav(3),yav(3),t(3) -c dimension a(3,3),b(3),c(3),d(3) -c do 10 j=1,3 -c d(j)=c(j) -c do 10 i=1,3 -c 10 d(j)=d(j)+a(j,i)*b(i) - do 10 j=1,3 - t(j)=yav(j) - do 10 i=1,3 - 10 t(j)=t(j)+b(j,i)*xav(i) - return - end - double precision function det (a,b,c) - implicit real*8 (a-h,o-z) - dimension a(3),b(3),c(3) - det=a(1)*(b(2)*c(3)-b(3)*c(2))+a(2)*(b(3)*c(1)-b(1)*c(3)) - 1 +a(3)*(b(1)*c(2)-b(2)*c(1)) - return - end - subroutine mmmul(a,b,c) - implicit real*8 (a-h,o-z) - dimension a(3,3),b(3,3),c(3,3) - do 10 i=1,3 - do 10 j=1,3 - c(i,j)=0.0d0 - do 10 k=1,3 - 10 c(i,j)=c(i,j)+a(i,k)*b(k,j) - return - end - subroutine matvec(uvec,tmat,pvec,nback) - implicit real*8 (a-h,o-z) - real*8 tmat(3,3),uvec(3,nback), pvec(3,nback) -c - do 2 j=1,nback - do 1 i=1,3 - uvec(i,j) = 0.0d0 - do 1 k=1,3 - 1 uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j) - 2 continue - return - end