X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=inline;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fdistfit.f;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fdistfit.f;h=80e8fe4736aa8ca53c5a1fb302a49f2c9a010b59;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/distfit.f b/source/unres/src_MD-M-newcorr/distfit.f new file mode 100644 index 0000000..80e8fe4 --- /dev/null +++ b/source/unres/src_MD-M-newcorr/distfit.f @@ -0,0 +1,207 @@ + subroutine distfit(debug,maxit) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.DISTFIT' + DIMENSION X(MAXRES),DIAGH(MAXRES),phiold(maxres) + logical debug,sing + +cinput------------------------------------ +c NX=NRES-3 +c NY=((NRES-4)*(NRES-5))/2 +cinput------------------------------------ +ctest MAXIT=20 + TOL=0.5 + MAXMAR=10 + RL=100.0 + + CALL TRANSFER(NRES,phi,phiold) + + F0=RDIF() + +cd WRITE (IOUT,*) 'DISTFIT: F0=',F0 + + + DO IT=1,MAXIT + CALL RDERIV + CALL HEVAL + + DO I=1,NX + DIAGH(I)=H(I,I) + ENDDO + RL=RL*0.1 + + DO IMAR=1,MAXMAR + DO I=1,NX + H(I,I)=DIAGH(I)+RL + ENDDO + CALL TRANSFER(NX,XX,X) + CALL BANACH(NX,MAXRES,H,X,sing) + AIN=0.0 + DO I=1,NX + AIN=AIN+DABS(X(I)) + ENDDO + IF (AIN.LT.0.1*TOL .AND. RL.LT.1.0E-4) THEN + if (debug) then + WRITE (IOUT,*) 'DISTFIT: CONVERGENCE HAS BEEN ACHIEVED' + WRITE (IOUT,*) 'IT=',it,'F=',F0 + endif + RETURN + ENDIF + DO I=4,NRES + phi(I)=phiold(I)+mask(i)*X(I-3) +c print *,X(I-3) + ENDDO + + F1=RDIF() +cd WRITE (IOUT,*) 'IMAR=',IMAR,' RL=',RL,' F1=',F1 + IF (F1.LT.F0) THEN + CALL TRANSFER(NRES,phi,phiold) + F0=F1 + GOTO 1 + ELSE IF (DABS(F1-F0).LT.1.0E-5) THEN + if (debug) then + WRITE (IOUT,*) 'DISTFIT: CANNOT IMPROVE DISTANCE FIT' + WRITE (IOUT,*) 'IT=',it,'F=',F1 + endif + RETURN + ENDIF + RL=RL*10.0 + ENDDO + WRITE (IOUT,*) 'DISTFIT: MARQUARDT PROCEDURE HAS FAILED' + WRITE (IOUT,*) 'IT=',it,'F=',F0 + CALL TRANSFER(NRES,phiold,phi) + RETURN + 1 continue +cd write (iout,*) "it",it," imar",imar," f0",f0 + enddo + WRITE (IOUT,*) 'DISTFIT: FINAL F=',F0,'after MAXIT=',maxit + return + END + + double precision FUNCTION RDIF() + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.DISTFIT' + +c print *,'in rdif' + + suma=0.0 + ind=0 + call chainbuild + do i=1,nres-3 + do j=i+3,nres + ind=ind+1 + if (w(ind).ne.0.0) then + DIJ=DIST(i,j) + suma=suma+w(ind)*(DIJ-d0(ind))*(DIJ-d0(ind)) + DD(ind)=DIJ +c print '(2i3,i4,4f12.2)',i,j,ind,dij,d0(ind),w(ind),suma + endif + enddo + enddo + + RDIF=suma + RETURN + END + + SUBROUTINE RDERIV + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.DISTFIT' + include 'COMMON.GEO' + DIMENSION E12(3),R13(3),R24(3),PRODU(3) + + DO I=1,NY + DO J=1,NX + DRDG(I,J)=0.0 + ENDDO + ENDDO + DO I=1,NX + I1=I+1 + I2=I+2 + CALL VEC(I1,I2,E12) + DO J=1,I + DO K=1,3 + R13(K)=C(K,J)-C(K,I1) + ENDDO + DO K=I2+1,NRES + DO L=1,3 + R24(L)=C(L,K)-C(L,I2) + ENDDO + IND=((J-1)*(2*NRES-J-6))/2+K-3 + PRODU(1)=R13(2)*R24(3)-R13(3)*R24(2) + PRODU(2)=R13(3)*R24(1)-R13(1)*R24(3) + PRODU(3)=R13(1)*R24(2)-R13(2)*R24(1) + DRDG(IND,I)=SCALAR(E12,PRODU)/DIST(J,K) + ENDDO + ENDDO + ENDDO + RETURN + END + + SUBROUTINE HEVAL + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.DISTFIT' + + DO I=1,NX + XI=0.0 + HII=0.0 + DO K=1,NY + BKI=DRDG(K,I) + BKIWK=w(K)*BKI + XI=XI+BKIWK*(D0(K)-DD(K)) + HII=HII+BKI*BKIWK + ENDDO + H(I,I)=HII + XX(I)=XI + DO J=I+1,NX + HIJ=0.0 + DO K=1,NY + HIJ=HIJ+DRDG(K,I)*DRDG(K,J)*w(K) + ENDDO + H(I,J)=HIJ + H(J,I)=HIJ + ENDDO + ENDDO + RETURN + END + + + SUBROUTINE VEC(I,J,U) +* +* Find the unit vector from atom (I) to atom (J). Store in U. +* + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + DIMENSION U(3) + + ANORM=0.0 + DO K=1,3 + UK=C(K,J)-C(K,I) + ANORM=ANORM+UK*UK + U(K)=UK + ENDDO + ANORM=SQRT(ANORM) + DO K=1,3 + U(K)=U(K)/ANORM + ENDDO + RETURN + END + + SUBROUTINE TRANSFER(N,X1,X2) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + DIMENSION X1(N),X2(N) + DO 1 I=1,N + 1 X2(I)=X1(I) + RETURN + END +