Fixed the homology gradient in finegrain mode
[unres.git] / source / wham / src / molread_zs.F
1       subroutine molread(*)
2 C
3 C Read molecular data.
4 C
5       implicit real*8 (a-h,o-z)
6       include 'DIMENSIONS'
7       include 'DIMENSIONS.ZSCOPT'
8       include 'COMMON.IOUNITS'
9       include 'COMMON.GEO'
10       include 'COMMON.VAR'
11 c     include 'include_unres/COMMON.VAR'
12       include 'COMMON.INTERACT'
13       include 'COMMON.LOCAL'
14       include 'COMMON.NAMES'
15       include 'COMMON.CHAIN'
16       include 'COMMON.FFIELD'
17       include 'COMMON.SBRIDGE'
18       include 'COMMON.TORCNSTR'
19       include 'COMMON.CONTROL'
20       character*4 sequence(maxres)
21       integer rescode
22       double precision x(maxvar)
23       character*320 controlcard,ucase
24       dimension itype_pdb(maxres)
25       logical seq_comp
26       call card_concat(controlcard,.true.)
27       call reada(controlcard,'SCAL14',scal14,0.4d0)
28       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
29       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
30       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
31 C     Bartek
32       call reada(controlcard,'WDFAD',wdfa_dist,0.0d0)
33       call reada(controlcard,'WDFAT',wdfa_tor,0.0d0)
34       call reada(controlcard,'WDFAN',wdfa_nei,0.0d0)
35       call reada(controlcard,'WDFAB',wdfa_beta,0.0d0)
36       write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
37      &  " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
38       r0_corr=cutoff_corr-delt_corr
39       call readi(controlcard,"NRES",nres,0)
40       iscode=index(controlcard,"ONE_LETTER")
41       if (nres.le.0) then
42         write (iout,*) "Error: no residues in molecule"
43         return1
44       endif
45       if (nres.gt.maxres) then
46         write (iout,*) "Error: too many residues",nres,maxres
47       endif
48       write(iout,*) 'nres=',nres
49 C Read sequence of the protein
50       if (iscode.gt.0) then
51         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
52       else
53         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
54       endif
55 C Convert sequence to numeric code
56       do i=1,nres
57         itype(i)=rescode(i,sequence(i),iscode)
58       enddo
59       write (iout,*) "Numeric code:"
60       write (iout,'(20i4)') (itype(i),i=1,nres)
61       do i=1,nres-1
62 #ifdef PROCOR
63         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
64 #else
65         if (itype(i).eq.21) then
66 #endif
67           itel(i)=0
68 #ifdef PROCOR
69         else if (itype(i+1).ne.20) then
70 #else
71         else if (itype(i).ne.20) then
72 #endif
73           itel(i)=1
74         else
75           itel(i)=2
76         endif
77       enddo
78       call read_bridge
79
80       if (with_dihed_constr) then
81
82       read (inp,*) ndih_constr
83       if (ndih_constr.gt.0) then
84         read (inp,*) ftors
85         write (iout,*) 'FTORS',ftors
86         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
87         write (iout,*)
88      &   'There are',ndih_constr,' constraints on phi angles.'
89         do i=1,ndih_constr
90           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
91         enddo
92         do i=1,ndih_constr
93           phi0(i)=deg2rad*phi0(i)
94           drange(i)=deg2rad*drange(i)
95         enddo
96       endif
97
98       endif
99
100       nnt=1
101       nct=nres
102       if (itype(1).eq.21) nnt=2
103       if (itype(nres).eq.21) nct=nct-1
104       write(iout,*) 'NNT=',NNT,' NCT=',NCT
105
106 C     Juyong:READ init_vars
107 C     Initialize variables!
108 C     Juyong:READ read_info
109 C     READ fragment information!!
110 C     both routines should be in dfa.F file!!
111
112       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
113      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
114        write (iout,*) "Calling init_dfa_vars"
115        call flush(iout)
116        call init_dfa_vars
117        write (iout,*) 'init_dfa_vars finished!'
118        call flush(iout)
119        call read_dfa_info
120        write (iout,*) 'read_dfa_info finished!'
121        call flush(iout)
122       endif
123
124 c Read distance restraints
125       if (constr_dist.gt.0) then
126         if (refstr) call read_ref_structure(*11)
127         call read_dist_constr
128         call hpb_partition
129       endif
130
131       if (constr_homology.gt.0) then
132 c       write (iout,*) "About to call read_constr_homology"
133 c       call flush(iout)
134         call read_constr_homology
135 c       write (iout,*) "Exit read_constr_homology"
136 c       call flush(iout)
137         if (indpdb.gt.0 .or. pdbref) then
138           do i=1,2*nres
139             do j=1,3
140               c(j,i)=crefjlee(j,i)
141               cref(j,i)=crefjlee(j,i)
142             enddo
143           enddo
144         endif
145 #ifdef DEBUG
146         write (iout,*) "Array C"
147         do i=1,nres
148           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
149      &      (c(j,i+nres),j=1,3)
150         enddo
151         write (iout,*) "Array Cref"
152         do i=1,nres
153           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
154      &      (cref(j,i+nres),j=1,3)
155         enddo
156 #endif
157 #ifdef DEBUG
158        call int_from_cart1(.false.)
159        call sc_loc_geom(.false.)
160        do i=1,nres
161          thetaref(i)=theta(i)
162          phiref(i)=phi(i)
163          write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
164        enddo
165        do i=1,nres-1
166          do j=1,3
167            dc(j,i)=c(j,i+1)-c(j,i)
168            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
169          enddo
170        enddo
171        do i=2,nres-1
172          do j=1,3
173            dc(j,i+nres)=c(j,i+nres)-c(j,i)
174            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
175          enddo
176        enddo
177 #endif
178       else
179         homol_nset=0
180       endif
181
182
183       call setup_var
184       call init_int_table
185       if (ns.gt.0) then
186         write (iout,'(/a,i3,a)') 'The chain contains',ns,
187      &  ' disulfide-bridging cysteines.'
188         write (iout,'(20i4)') (iss(i),i=1,ns)
189         write (iout,'(/a/)') 'Pre-formed links are:' 
190         do i=1,nss
191           i1=ihpb(i)-nres
192           i2=jhpb(i)-nres
193           it1=itype(i1)
194           it2=itype(i2)
195          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
196      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
197      &    dhpb(i),ebr,forcon(i)
198         enddo
199       endif
200       write (iout,'(a)')
201       return
202    11 stop "Error reading reference structure"
203       end
204 c-----------------------------------------------------------------------------
205       logical function seq_comp(itypea,itypeb,length)
206       implicit none
207       integer length,itypea(length),itypeb(length)
208       integer i
209       do i=1,length
210         if (itypea(i).ne.itypeb(i)) then
211           seq_comp=.false.
212           return
213         endif
214       enddo
215       seq_comp=.true.
216       return
217       end
218 c-----------------------------------------------------------------------------
219       subroutine read_bridge
220 C Read information about disulfide bridges.
221       implicit real*8 (a-h,o-z)
222       include 'DIMENSIONS'
223       include 'DIMENSIONS.ZSCOPT'
224       include 'COMMON.IOUNITS'
225       include 'COMMON.GEO'
226       include 'COMMON.VAR'
227       include 'COMMON.INTERACT'
228       include 'COMMON.NAMES'
229       include 'COMMON.CHAIN'
230       include 'COMMON.FFIELD'
231       include 'COMMON.SBRIDGE'
232 C Read bridging residues.
233       read (inp,*) ns,(iss(i),i=1,ns)
234       print *,'ns=',ns
235       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
236 C Check whether the specified bridging residues are cystines.
237       do i=1,ns
238         if (itype(iss(i)).ne.1) then
239           write (iout,'(2a,i3,a)') 
240      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
241      &   ' can form a disulfide bridge?!!!'
242           write (*,'(2a,i3,a)') 
243      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
244      &   ' can form a disulfide bridge?!!!'
245          stop
246         endif
247       enddo
248 C Read preformed bridges.
249       if (ns.gt.0) then
250       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
251       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
252       if (nss.gt.0) then
253         nhpb=nss
254 C Check if the residues involved in bridges are in the specified list of
255 C bridging residues.
256         do i=1,nss
257           do j=1,i-1
258             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
259      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
260               write (iout,'(a,i3,a)') 'Disulfide pair',i,
261      &      ' contains residues present in other pairs.'
262               write (*,'(a,i3,a)') 'Disulfide pair',i,
263      &      ' contains residues present in other pairs.'
264               stop 
265             endif
266           enddo
267           do j=1,ns
268             if (ihpb(i).eq.iss(j)) goto 10
269           enddo
270           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
271    10     continue
272           do j=1,ns
273             if (jhpb(i).eq.iss(j)) goto 20
274           enddo
275           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
276    20     continue
277           dhpb(i)=dbr
278           forcon(i)=fbr
279         enddo
280         do i=1,nss
281           ihpb(i)=ihpb(i)+nres
282           jhpb(i)=jhpb(i)+nres
283         enddo
284       endif
285       endif
286       if (ns.gt.0.and.dyn_ss) then
287 C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of
288 C the bond
289           do i=nss+1,nhpb
290 C /06/28/2013 Adasko: nss number of full SS bonds
291             ihpb(i-nss)=ihpb(i)
292             jhpb(i-nss)=jhpb(i)
293             forcon(i-nss)=forcon(i)
294             dhpb(i-nss)=dhpb(i)
295           enddo
296           nhpb=nhpb-nss
297           nss=0
298           call hpb_partition
299           do i=1,ns
300             dyn_ss_mask(iss(i))=.true.
301 C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond
302 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
303           enddo
304       endif
305       return
306       end
307 c------------------------------------------------------------------------------
308       subroutine read_angles(kanal,iscor,energ,iprot,*)
309       implicit real*8 (a-h,o-z)
310       include 'DIMENSIONS'
311       include 'DIMENSIONS.ZSCOPT'
312       include 'COMMON.INTERACT'
313       include 'COMMON.SBRIDGE'
314       include 'COMMON.GEO'
315       include 'COMMON.VAR'
316       include 'COMMON.CHAIN'
317       include 'COMMON.IOUNITS'
318       character*80 lineh
319       read(kanal,'(a80)',end=10,err=10) lineh
320       read(lineh(:5),*,err=8) ic
321       read(lineh(6:),*,err=8) energ
322       goto 9
323     8 ic=1
324       print *,'error, assuming e=1d10',lineh
325       energ=1d10
326       nss=0
327     9 continue
328       read(lineh(18:),*,end=10,err=10) nss
329       IF (NSS.LT.9) THEN
330         read (lineh(20:),*,end=10,err=10)
331      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
332       ELSE
333         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
334         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
335      &    I=9,NSS),iscor
336       ENDIF
337 c      print *,"energy",energ," iscor",iscor
338       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
339       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
340       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
341       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
342       do i=1,nres
343         theta(i)=deg2rad*theta(i)
344         phi(i)=deg2rad*phi(i)
345         alph(i)=deg2rad*alph(i)
346         omeg(i)=deg2rad*omeg(i)
347       enddo
348       return
349    10 return1
350       end
351 c-------------------------------------------------------------------------------
352       subroutine read_dist_constr
353       implicit real*8 (a-h,o-z)
354       include 'DIMENSIONS'
355       include 'DIMENSIONS.ZSCOPT'
356       include 'COMMON.CONTROL'
357       include 'COMMON.CHAIN'
358       include 'COMMON.IOUNITS'
359       include 'COMMON.SBRIDGE'
360       integer ifrag_(2,100),ipair_(2,100)
361       double precision wfrag_(100),wpair_(100)
362       character*500 controlcard
363 c      write (iout,*) "Calling read_dist_constr"
364 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
365 c      call flush(iout)
366       call card_concat(controlcard,.true.)
367       call readi(controlcard,"NFRAG",nfrag_,0)
368       call readi(controlcard,"NPAIR",npair_,0)
369       call readi(controlcard,"NDIST",ndist_,0)
370       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
371       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
372       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
373       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
374       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
375       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
376       write (iout,*) "IFRAG"
377       do i=1,nfrag_
378         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
379       enddo
380       write (iout,*) "IPAIR"
381       do i=1,npair_
382         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
383       enddo
384       call flush(iout)
385       if (.not.refstr .and. nfrag_.gt.0) then
386         write (iout,*) 
387      &  "ERROR: no reference structure to compute distance restraints"
388         write (iout,*)
389      &  "Restraints must be specified explicitly (NDIST=number)"
390         stop 
391       endif
392       if (nfrag_.lt.2 .and. npair_.gt.0) then 
393         write (iout,*) "ERROR: Less than 2 fragments specified",
394      &   " but distance restraints between pairs requested"
395         stop 
396       endif 
397       call flush(iout)
398       do i=1,nfrag_
399         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
400         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
401      &    ifrag_(2,i)=nstart_sup+nsup-1
402 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
403         call flush(iout)
404         if (wfrag_(i).gt.0.0d0) then
405         do j=ifrag_(1,i),ifrag_(2,i)-1
406           do k=j+1,ifrag_(2,i)
407             write (iout,*) "j",j," k",k
408             ddjk=dist(j,k)
409             if (constr_dist.eq.1) then
410             nhpb=nhpb+1
411             ihpb(nhpb)=j
412             jhpb(nhpb)=k
413               dhpb(nhpb)=ddjk
414             forcon(nhpb)=wfrag_(i) 
415             else if (constr_dist.eq.2) then
416               if (ddjk.le.dist_cut) then
417                 nhpb=nhpb+1
418                 ihpb(nhpb)=j
419                 jhpb(nhpb)=k
420                 dhpb(nhpb)=ddjk
421                 forcon(nhpb)=wfrag_(i) 
422               endif
423             else
424               nhpb=nhpb+1
425               ihpb(nhpb)=j
426               jhpb(nhpb)=k
427               dhpb(nhpb)=ddjk
428               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
429             endif
430             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
431      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
432           enddo
433         enddo
434         endif
435       enddo
436       do i=1,npair_
437         if (wpair_(i).gt.0.0d0) then
438         ii = ipair_(1,i)
439         jj = ipair_(2,i)
440         if (ii.gt.jj) then
441           itemp=ii
442           ii=jj
443           jj=itemp
444         endif
445         do j=ifrag_(1,ii),ifrag_(2,ii)
446           do k=ifrag_(1,jj),ifrag_(2,jj)
447             nhpb=nhpb+1
448             ihpb(nhpb)=j
449             jhpb(nhpb)=k
450             forcon(nhpb)=wpair_(i)
451             dhpb(nhpb)=dist(j,k)
452             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
453      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
454           enddo
455         enddo
456         endif
457       enddo 
458       do i=1,ndist_
459         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
460      &     ibecarb(i),forcon(nhpb+1)
461         if (forcon(nhpb+1).gt.0.0d0) then
462           nhpb=nhpb+1
463           if (ibecarb(i).gt.0) then
464             ihpb(i)=ihpb(i)+nres
465             jhpb(i)=jhpb(i)+nres
466           endif
467           if (dhpb(nhpb).eq.0.0d0) 
468      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
469         endif
470       enddo
471       do i=1,nhpb
472           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
473      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
474       enddo
475       call flush(iout)
476       return
477       end
478
479
480
481 c====-------------------------------------------------------------------
482       subroutine read_constr_homology
483
484       include 'DIMENSIONS'
485       include 'DIMENSIONS.ZSCOPT'
486 #ifdef MPI
487       include 'mpif.h'
488 #endif
489       include 'COMMON.SETUP'
490       include 'COMMON.CONTROL'
491       include 'COMMON.CHAIN'
492       include 'COMMON.IOUNITS'
493       include 'COMMON.GEO'
494       include 'COMMON.INTERACT'
495       include 'COMMON.HOMRESTR'
496 c
497 c For new homol impl
498 c
499       include 'COMMON.VAR'
500 c     include 'include_unres/COMMON.VAR'
501 c
502
503 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
504 c    &                 dist_cut
505 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
506 c    &    sigma_odl_temp(maxres,maxres,max_template)
507       character*2 kic2
508       character*24 model_ki_dist, model_ki_angle
509       character*500 controlcard
510       integer ki, i, j, k, l
511       logical lprn /.true./
512 c
513 c     FP - Nov. 2014 Temporary specifications for new vars
514 c
515       double precision rescore_tmp,x12,y12,z12
516       double precision, dimension (max_template,maxres) :: rescore
517       character*24 tpl_k_rescore
518 c -----------------------------------------------------------------
519 c Reading multiple PDB ref structures and calculation of retraints
520 c not using pre-computed ones stored in files model_ki_{dist,angle}
521 c FP (Nov., 2014)
522 c -----------------------------------------------------------------
523 c
524 c
525 c Alternative: reading from input
526       call card_concat(controlcard,.true.)
527       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
528       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
529       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
530       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
531       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
532  
533       call readi(controlcard,"HOMOL_NSET",homol_nset,1)       
534       if (homol_nset.gt.1)then
535          call card_concat(controlcard,.true.)
536          read(controlcard,*) (waga_homology(i),i=1,homol_nset) 
537          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
538           write(iout,*) "iset homology_weight "
539 c         do i=1,homol_nset
540 c          write(iout,*) i,waga_homology(i)
541 c         enddo
542          endif
543          iset=mod(kolor,homol_nset)+1
544       else
545        iset=1
546        waga_homology(1)=1.0
547       endif
548 c     write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
549
550 cd      write (iout,*) "nnt",nnt," nct",nct
551 cd      call flush(iout)
552
553
554       lim_odl=0
555       lim_dih=0
556 c
557 c  New
558 c
559       lim_theta=0
560       lim_xx=0
561 c
562 c  Reading HM global scores (prob not required)
563 c
564 c     open (4,file="HMscore")
565 c     do k=1,constr_homology
566 c       read (4,*,end=521) hmscore_tmp
567 c       hmscore(k)=hmscore_tmp ! Another transformation can be used 
568 c       write(*,*) "Model", k, ":", hmscore(k)
569 c     enddo
570 c521  continue
571
572 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
573
574       do k=1,constr_homology
575
576         read(inp,'(a)') pdbfile
577 c  Next stament causes error upon compilation (?)
578 c       if(me.eq.king.or. .not. out1file)
579 c         write (iout,'(2a)') 'PDB data will be read from file ',
580 c    &   pdbfile(:ilen(pdbfile))
581         open(ipdbin,file=pdbfile,status='old',err=33)
582         goto 34
583   33    write (iout,'(a)') 'Error opening PDB file.'
584         stop
585   34    continue
586 c        print *,'Begin reading pdb data'
587 c
588 c Files containing res sim or local scores (former containing sigmas)
589 c
590
591         write(kic2,'(bz,i2.2)') k
592
593         tpl_k_rescore="template"//kic2//".sco"
594 c       tpl_k_sigma_odl="template"//kic2//".sigma_odl"
595 c       tpl_k_sigma_dih="template"//kic2//".sigma_dih"
596 c       tpl_k_sigma_theta="template"//kic2//".sigma_theta"
597 c       tpl_k_sigma_d="template"//kic2//".sigma_d"
598
599         unres_pdb=.false.
600         call readpdb
601         do i=1,2*nres
602           do j=1,3
603             crefjlee(j,i)=c(j,i)
604           enddo
605         enddo
606 #ifdef DEBUG
607         do i=1,nres
608           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
609      &      (crefjlee(j,i+nres),j=1,3)
610         enddo
611 #endif
612         write (iout,*) "read_constr_homology: after reading pdb file"
613         call flush(iout)
614
615 c
616 c     Distance restraints
617 c
618 c          ... --> odl(k,ii)
619 C Copy the coordinates from reference coordinates (?)
620         do i=1,2*nres
621           do j=1,3
622             c(j,i)=cref(j,i)
623 c           write (iout,*) "c(",j,i,") =",c(j,i)
624           enddo
625         enddo
626 c
627 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
628 c
629 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
630           open (ientin,file=tpl_k_rescore,status='old')
631           do irec=1,maxdim ! loop for reading res sim 
632             if (irec.eq.1) then
633                rescore(k,irec)=0.0d0
634                goto 1301 
635             endif
636             read (ientin,*,end=1401) rescore_tmp
637 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
638             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
639 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
640  1301     continue
641           enddo  
642  1401   continue
643           close (ientin)        
644 c         open (ientin,file=tpl_k_sigma_odl,status='old')
645 c         do irec=1,maxdim ! loop for reading sigma_odl
646 c            read (ientin,*,end=1401) i, j, 
647 c    &                                sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
648 c            sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
649 c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k) 
650 c         enddo
651 c 1401   continue
652 c         close (ientin)
653         if (waga_dist.ne.0.0d0) then
654           ii=0
655           do i = nnt,nct-2 ! right? without parallel.
656             do j=i+2,nct ! right?
657 c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology 
658 c           do j=i+2,nres ! ibid
659 c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
660 c           do j=i+2,nct ! ibid
661               ii=ii+1
662 c             write (iout,*) "k",k
663 c             write (iout,*) "i",i," j",j," constr_homology",
664 c    &                       constr_homology
665               ires_homo(ii)=i
666               jres_homo(ii)=j
667 c
668 c Attempt to replace dist(i,j) by its definition in ...
669 c
670               x12=c(1,i)-c(1,j)
671               y12=c(2,i)-c(2,j)
672               z12=c(3,i)-c(3,j)
673               distal=dsqrt(x12*x12+y12*y12+z12*z12)
674               odl(k,ii)=distal
675 c
676 c             odl(k,ii)=dist(i,j)
677 c             write (iout,*) "dist(",i,j,") =",dist(i,j)
678 c             write (iout,*) "distal = ",distal
679 c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
680 c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
681 c    &                      "rescore(",k,j,") =",rescore(k,j)
682 c
683 c  Calculation of sigma from res sim
684 c
685 c             if (odl(k,ii).le.6.0d0) then
686 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
687 c    Other functional forms possible depending on odl(k,ii), eg.
688 c
689             if (odl(k,ii).le.dist_cut) then
690               sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
691 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
692             else
693 #ifdef OLDSIGMA
694               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
695      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
696 #else
697               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
698      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
699 #endif
700 c   Following expr replaced by a positive exp argument
701 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
702 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
703
704 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
705 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
706             endif
707 c
708               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
709 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
710 c
711 c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
712 c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
713             enddo
714 c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
715           enddo
716 c         lim_odl=ii
717 c         if (constr_homology.gt.0) call homology_partition
718         endif
719 c
720 c     Theta, dihedral and SC retraints
721 c
722         if (waga_angle.gt.0.0d0) then
723 c         open (ientin,file=tpl_k_sigma_dih,status='old')
724 c         do irec=1,maxres-3 ! loop for reading sigma_dih
725 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
726 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
727 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
728 c    &                            sigma_dih(k,i+nnt-1)
729 c         enddo
730 c1402   continue
731 c         close (ientin)
732           do i = nnt+3,nct ! right? without parallel.
733 c         do i=1,nres ! alternative for bounds acc to readpdb?
734 c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
735 c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
736             dih(k,i)=phiref(i) ! right?
737 c           read (ientin,*) sigma_dih(k,i) ! original variant
738 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
739 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
740 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
741 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
742 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
743
744             sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
745      &                     rescore(k,i-2)+rescore(k,i-3)  !  right expression ?
746 c
747 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
748 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
749 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
750 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
751 c   Instead of res sim other local measure of b/b str reliability possible
752             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
753 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
754             if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
755 c           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
756           enddo
757         endif
758
759         if (waga_theta.gt.0.0d0) then
760 c         open (ientin,file=tpl_k_sigma_theta,status='old')
761 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
762 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
763 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
764 c    &                              sigma_theta(k,i+nnt-1)
765 c         enddo
766 c1403   continue
767 c         close (ientin)
768
769           do i = nnt+2,nct ! right? without parallel.
770 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
771 c         do i=ithet_start,ithet_end ! with FG parallel.
772              thetatpl(k,i)=thetaref(i)
773 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
774 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
775 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
776 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
777 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
778              sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
779      &                        rescore(k,i-2) !  right expression ?
780              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
781
782 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
783 c                             rescore(k,i-2) !  right expression ?
784 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
785              if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
786           enddo
787         endif
788
789         if (waga_d.gt.0.0d0) then
790 c       open (ientin,file=tpl_k_sigma_d,status='old')
791 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
792 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
793 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
794 c    &                          sigma_d(k,i+nnt-1)
795 c         enddo
796 c1404   continue
797           close (ientin)
798
799           do i = nnt,nct ! right? without parallel.
800 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
801 c         do i=loc_start,loc_end ! with FG parallel.
802              if (itype(i).eq.10) goto 1 ! right?
803                xxtpl(k,i)=xxref(i)
804                yytpl(k,i)=yyref(i)
805                zztpl(k,i)=zzref(i)
806 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
807 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
808 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
809 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
810                sigma_d(k,i)=rescore(k,i) !  right expression ?
811                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
812
813 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
814 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
815 c              read (ientin,*) sigma_d(k,i) ! 1st variant
816                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
817     1     continue
818           enddo
819         endif
820         close(ientin)
821       enddo
822       if (waga_dist.ne.0.0d0) lim_odl=ii
823       if (constr_homology.gt.0) call homology_partition
824       if (constr_homology.gt.0) call init_int_table
825 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
826 cd     & "lim_xx=",lim_xx
827 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
828 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
829 c
830 c Print restraints
831 c
832       if (.not.lprn) return
833 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
834       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
835        write (iout,*) "Distance restraints from templates"
836        do ii=1,lim_odl
837        write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
838      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
839        enddo
840        write (iout,*) "Dihedral angle restraints from templates"
841        do i=nnt+3,lim_dih
842         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
843      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
844        enddo
845        write (iout,*) "Virtual-bond angle restraints from templates"
846        do i=nnt+2,lim_theta
847         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
848      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
849        enddo
850        write (iout,*) "SC restraints from templates"
851        do i=nnt,lim_xx
852         write(iout,'(i5,10(4f8.2,4x))') i,
853      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
854      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
855        enddo
856       endif
857 c -----------------------------------------------------------------
858       return
859       end
860 c----------------------------------------------------------------------