80e8fe4736aa8ca53c5a1fb302a49f2c9a010b59
[unres.git] / source / unres / src_MD / 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