X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Ffitsq.f;fp=source%2Fwham%2Fsrc-NEWSC%2Ffitsq.f;h=17d92eef5f8483be8ede0bbb06726d69ae0521c3;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/wham/src-NEWSC/fitsq.f b/source/wham/src-NEWSC/fitsq.f new file mode 100755 index 0000000..17d92ee --- /dev/null +++ b/source/wham/src-NEWSC/fitsq.f @@ -0,0 +1,352 @@ + 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 + 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 +crc 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 + 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 + eta = z00100000 + nit = 0 + small=25.0*10.e-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 + 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 + if (np.gt.npq) go to 230 + n0=0 + do 220 n=np,npq + nn=n+np-1 + 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 + 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 + 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 + 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 + 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 + 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)) +* write (1,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