f359cac805901b77164844d05903abe3ca850ac7
[unres.git] / source / wham / src-M / read_dist_constr.F
1       subroutine read_dist_constr
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'DIMENSIONS.FREE'
5       include 'COMMON.CONTROL'
6       include 'COMMON.CHAIN'
7       include 'COMMON.IOUNITS'
8       include 'COMMON.SBRIDGE'
9       integer ifrag_(2,100),ipair_(2,100)
10       double precision wfrag_(100),wpair_(100)
11       character*500 controlcard
12       logical normalize,next
13       integer restr_type
14       double precision xlink(4,0:4) /
15 c           a          b       c     sigma
16      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
17      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
18      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
19      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
20      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
21       write (iout,*) "Calling read_dist_constr"
22 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
23 c      call flush(iout)
24       next=.true.
25
26       DO WHILE (next)
27
28       call card_concat(controlcard,.true.)
29       next = index(controlcard,"NEXT").gt.0
30       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
31       write (iout,*) "restr_type",restr_type
32       call readi(controlcard,"NFRAG",nfrag_,0)
33       call readi(controlcard,"NFRAG",nfrag_,0)
34       call readi(controlcard,"NPAIR",npair_,0)
35       call readi(controlcard,"NDIST",ndist_,0)
36       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
37       if (restr_type.eq.10) 
38      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
39       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
40       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
41       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
42       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
43       normalize = index(controlcard,"NORMALIZE").gt.0
44       write (iout,*) "WBOLTZD",wboltzd
45 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
46 c      write (iout,*) "IFRAG"
47 c      do i=1,nfrag_
48 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
49 c      enddo
50 c      write (iout,*) "IPAIR"
51 c      do i=1,npair_
52 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
53 c      enddo
54       if (nfrag_.gt.0) then
55         nres0=nres
56         read(inp,'(a)') pdbfile
57         write (iout,*) 
58      & "Distance restraints will be constructed from structure ",pdbfile
59         open(ipdbin,file=pdbfile,status='old',err=11)
60         call readpdb(.true.)
61         nres=nres0
62         close(ipdbin)
63       endif
64       do i=1,nfrag_
65         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
66         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
67      &    ifrag_(2,i)=nstart_sup+nsup-1
68 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
69 c        call flush(iout)
70         if (wfrag_(i).eq.0.0d0) cycle
71         do j=ifrag_(1,i),ifrag_(2,i)-1
72           do k=j+1,ifrag_(2,i)
73 c            write (iout,*) "j",j," k",k
74             ddjk=dist(j,k)
75             if (restr_type.eq.1) then
76               nhpb=nhpb+1
77               irestr_type(nhpb)=1
78               ihpb(nhpb)=j
79               jhpb(nhpb)=k
80               dhpb(nhpb)=ddjk
81               forcon(nhpb)=wfrag_(i) 
82             else if (constr_dist.eq.2) then
83               if (ddjk.le.dist_cut) then
84                 nhpb=nhpb+1
85                 irestr_type(nhpb)=1
86                 ihpb(nhpb)=j
87                 jhpb(nhpb)=k
88                 dhpb(nhpb)=ddjk
89                 forcon(nhpb)=wfrag_(i) 
90               endif
91             else if (restr_type.eq.3) then
92               nhpb=nhpb+1
93               irestr_type(nhpb)=1
94               ihpb(nhpb)=j
95               jhpb(nhpb)=k
96               dhpb(nhpb)=ddjk
97               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
98             endif
99             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
100      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
101           enddo
102         enddo
103       enddo
104       do i=1,npair_
105         if (wpair_(i).eq.0.0d0) cycle
106         ii = ipair_(1,i)
107         jj = ipair_(2,i)
108         if (ii.gt.jj) then
109           itemp=ii
110           ii=jj
111           jj=itemp
112         endif
113         do j=ifrag_(1,ii),ifrag_(2,ii)
114           do k=ifrag_(1,jj),ifrag_(2,jj)
115             if (restr_type.eq.1) then
116               nhpb=nhpb+1
117               irestr_type(nhpb)=1
118               ihpb(nhpb)=j
119               jhpb(nhpb)=k
120               dhpb(nhpb)=ddjk
121               forcon(nhpb)=wfrag_(i) 
122             else if (constr_dist.eq.2) then
123               if (ddjk.le.dist_cut) then
124                 nhpb=nhpb+1
125                 irestr_type(nhpb)=1
126                 ihpb(nhpb)=j
127                 jhpb(nhpb)=k
128                 dhpb(nhpb)=ddjk
129                 forcon(nhpb)=wfrag_(i) 
130               endif
131             else if (restr_type.eq.3) then
132               nhpb=nhpb+1
133               irestr_type(nhpb)=1
134               ihpb(nhpb)=j
135               jhpb(nhpb)=k
136               dhpb(nhpb)=ddjk
137               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
138             endif
139             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
140      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
141           enddo
142         enddo
143       enddo 
144
145 c      print *,ndist_
146       write (iout,*) "Distance restraints as read from input"
147       do i=1,ndist_
148         if (restr_type.eq.11) then
149           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
150      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
151 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
152           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
153           nhpb=nhpb+1
154           irestr_type(nhpb)=11
155           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
156      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
157      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
158           if (ibecarb(nhpb).gt.0) then
159             ihpb(nhpb)=ihpb(nhpb)+nres
160             jhpb(nhpb)=jhpb(nhpb)+nres
161           endif
162         else if (constr_dist.eq.10) then
163 c Cross-lonk Markov-like potential
164           call card_concat(controlcard,.true.)
165           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
166           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
167           ibecarb(nhpb+1)=0
168           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
169           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
170           if (index(controlcard,"ZL").gt.0) then
171             link_type=1
172           else if (index(controlcard,"ADH").gt.0) then
173             link_type=2
174           else if (index(controlcard,"PDH").gt.0) then
175             link_type=3
176           else if (index(controlcard,"DSS").gt.0) then
177             link_type=4
178           else
179             link_type=0
180           endif
181           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
182      &       xlink(1,link_type))
183           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
184      &       xlink(2,link_type))
185           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
186      &       xlink(3,link_type))
187           call reada(controlcard,"SIGMA",forcon(nhpb+1),
188      &       xlink(4,link_type))
189           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
190 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
191 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
192           if (forcon(nhpb+1).le.0.0d0 .or. 
193      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
194           nhpb=nhpb+1
195           irestr_type(nhpb)=10
196           if (ibecarb(nhpb).gt.0) then
197             ihpb(nhpb)=ihpb(nhpb)+nres
198             jhpb(nhpb)=jhpb(nhpb)+nres
199           endif
200           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
201      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
202      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
203      &     irestr_type(nhpb)
204         else
205 C        print *,"in else"
206           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
207      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
208           if (forcon(nhpb+1).gt.0.0d0) then
209           nhpb=nhpb+1
210           if (dhpb1(nhpb).eq.0.0d0) then
211             irestr_type(nhpb)=1
212           else
213             irestr_type(nhpb)=2
214           endif
215           if (ibecarb(nhpb).gt.0) then
216             ihpb(nhpb)=ihpb(nhpb)+nres
217             jhpb(nhpb)=jhpb(nhpb)+nres
218           endif
219           if (dhpb(nhpb).eq.0.0d0)
220      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
221         endif
222         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
223      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
224         endif
225 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
226 C        if (forcon(nhpb+1).gt.0.0d0) then
227 C          nhpb=nhpb+1
228 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
229       enddo
230
231       ENDDO ! next
232
233       fordepthmax=0.0d0
234       if (normalize) then
235         do i=nss+1,nhpb
236           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax) 
237      &      fordepthmax=fordepth(i)
238         enddo
239         do i=nss+1,nhpb
240           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
241         enddo
242       endif
243       if (nhpb.gt.nss)  then
244         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
245      &  "The following",nhpb-nss,
246      &  " distance restraints have been imposed:",
247      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
248      &  "  score"," type"
249         do i=nss+1,nhpb
250           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
251      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
252      &  irestr_type(i)
253         enddo
254       endif
255       write (iout,*) 
256       call hpb_partition
257       call flush(iout)
258       return
259    11 write (iout,*)"read_dist_restr: error reading reference structure"
260       stop
261       end