source/unres/src_MD-M/distfit.f

   1       subroutine distfit(debug,maxit)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.CHAIN'
   5       include 'COMMON.VAR'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.DISTFIT'
   8       DIMENSION X(MAXRES),DIAGH(MAXRES),phiold(maxres)
   9       logical debug,sing
  10
  11 cinput------------------------------------
  12 c       NX=NRES-3
  13 c       NY=((NRES-4)*(NRES-5))/2
  14 cinput------------------------------------
  15 ctest      MAXIT=20
  16       TOL=0.5
  17       MAXMAR=10
  18       RL=100.0
  19
  20       CALL TRANSFER(NRES,phi,phiold)
  21
  22       F0=RDIF()
  23
  24 cd      WRITE (IOUT,*) 'DISTFIT: F0=',F0
  25
  26
  27       DO IT=1,MAXIT
  28         CALL RDERIV
  29         CALL HEVAL
  30
  31         DO I=1,NX
  32           DIAGH(I)=H(I,I)
  33         ENDDO
  34         RL=RL*0.1
  35
  36         DO IMAR=1,MAXMAR
  37           DO I=1,NX
  38             H(I,I)=DIAGH(I)+RL
  39           ENDDO
  40           CALL TRANSFER(NX,XX,X)
  41           CALL BANACH(NX,MAXRES,H,X,sing)
  42           AIN=0.0
  43           DO I=1,NX
  44             AIN=AIN+DABS(X(I))
  45           ENDDO
  46           IF (AIN.LT.0.1*TOL .AND. RL.LT.1.0E-4) THEN
  47               if (debug) then
  48               WRITE (IOUT,*) 'DISTFIT: CONVERGENCE HAS BEEN ACHIEVED'
  49               WRITE (IOUT,*) 'IT=',it,'F=',F0
  50               endif
  51             RETURN
  52           ENDIF
  53           DO I=4,NRES
  54             phi(I)=phiold(I)+mask(i)*X(I-3)
  55 c            print *,X(I-3)
  56           ENDDO
  57
  58           F1=RDIF()
  59 cd          WRITE (IOUT,*) 'IMAR=',IMAR,' RL=',RL,' F1=',F1
  60           IF (F1.LT.F0) THEN
  61             CALL TRANSFER(NRES,phi,phiold)
  62             F0=F1
  63             GOTO 1
  64           ELSE IF (DABS(F1-F0).LT.1.0E-5) THEN
  65             if (debug) then
  66             WRITE (IOUT,*) 'DISTFIT: CANNOT IMPROVE DISTANCE FIT'
  67             WRITE (IOUT,*) 'IT=',it,'F=',F1
  68             endif
  69             RETURN
  70           ENDIF
  71           RL=RL*10.0
  72         ENDDO
  73         WRITE (IOUT,*) 'DISTFIT: MARQUARDT PROCEDURE HAS FAILED'
  74         WRITE (IOUT,*) 'IT=',it,'F=',F0
  75         CALL TRANSFER(NRES,phiold,phi)
  76         RETURN
  77     1   continue
  78 cd        write (iout,*) "it",it," imar",imar," f0",f0
  79       enddo
  80       WRITE (IOUT,*) 'DISTFIT: FINAL F=',F0,'after MAXIT=',maxit
  81       return
  82       END
  83
  84       double precision FUNCTION RDIF()
  85       implicit real*8 (a-h,o-z)
  86       include 'DIMENSIONS'
  87       include 'COMMON.CHAIN'
  88       include 'COMMON.DISTFIT'
  89
  90 c      print *,'in rdif'
  91
  92       suma=0.0
  93       ind=0
  94       call chainbuild
  95       do i=1,nres-3
  96         do j=i+3,nres
  97           ind=ind+1
  98           if (w(ind).ne.0.0) then
  99             DIJ=DIST(i,j)
 100             suma=suma+w(ind)*(DIJ-d0(ind))*(DIJ-d0(ind))
 101             DD(ind)=DIJ
 102 c            print '(2i3,i4,4f12.2)',i,j,ind,dij,d0(ind),w(ind),suma
 103           endif
 104         enddo
 105       enddo
 106
 107       RDIF=suma
 108       RETURN
 109       END
 110
 111       SUBROUTINE RDERIV
 112       implicit real*8 (a-h,o-z)
 113       include 'DIMENSIONS'
 114       include 'COMMON.CHAIN'
 115       include 'COMMON.DISTFIT'
 116       include 'COMMON.GEO'
 117       DIMENSION E12(3),R13(3),R24(3),PRODU(3)
 118
 119       DO I=1,NY
 120         DO J=1,NX
 121           DRDG(I,J)=0.0
 122         ENDDO
 123       ENDDO
 124       DO I=1,NX
 125         I1=I+1
 126         I2=I+2
 127         CALL VEC(I1,I2,E12)
 128         DO J=1,I
 129           DO K=1,3
 130             R13(K)=C(K,J)-C(K,I1)
 131           ENDDO
 132           DO K=I2+1,NRES
 133             DO L=1,3
 134               R24(L)=C(L,K)-C(L,I2)
 135             ENDDO
 136             IND=((J-1)*(2*NRES-J-6))/2+K-3
 137             PRODU(1)=R13(2)*R24(3)-R13(3)*R24(2)
 138             PRODU(2)=R13(3)*R24(1)-R13(1)*R24(3)
 139             PRODU(3)=R13(1)*R24(2)-R13(2)*R24(1)
 140             DRDG(IND,I)=SCALAR(E12,PRODU)/DIST(J,K)
 141           ENDDO
 142         ENDDO
 143       ENDDO
 144       RETURN
 145       END
 146
 147       SUBROUTINE HEVAL
 148       implicit real*8 (a-h,o-z)
 149       include 'DIMENSIONS'
 150       include 'COMMON.CHAIN'
 151       include 'COMMON.DISTFIT'
 152
 153       DO I=1,NX
 154         XI=0.0
 155         HII=0.0
 156         DO K=1,NY
 157           BKI=DRDG(K,I)
 158           BKIWK=w(K)*BKI
 159           XI=XI+BKIWK*(D0(K)-DD(K))
 160           HII=HII+BKI*BKIWK
 161         ENDDO
 162         H(I,I)=HII
 163         XX(I)=XI
 164         DO J=I+1,NX
 165           HIJ=0.0
 166           DO K=1,NY
 167             HIJ=HIJ+DRDG(K,I)*DRDG(K,J)*w(K)
 168           ENDDO
 169           H(I,J)=HIJ
 170           H(J,I)=HIJ
 171         ENDDO
 172       ENDDO
 173       RETURN
 174       END
 175
 176
 177       SUBROUTINE VEC(I,J,U)
 178 *
 179 *  Find the unit vector from atom (I) to atom (J). Store in U.
 180 *
 181       implicit real*8 (a-h,o-z)
 182       include 'DIMENSIONS'
 183       include 'COMMON.CHAIN'
 184       DIMENSION U(3)
 185
 186       ANORM=0.0
 187       DO K=1,3
 188         UK=C(K,J)-C(K,I)
 189         ANORM=ANORM+UK*UK
 190         U(K)=UK
 191       ENDDO
 192       ANORM=SQRT(ANORM)
 193       DO K=1,3
 194         U(K)=U(K)/ANORM
 195       ENDDO
 196       RETURN
 197       END
 198
 199       SUBROUTINE TRANSFER(N,X1,X2)
 200       implicit real*8 (a-h,o-z)
 201       include 'DIMENSIONS'
 202       DIMENSION X1(N),X2(N)
 203       DO 1 I=1,N
 204     1   X2(I)=X1(I)
 205       RETURN
 206       END
 207