correction to the last commit - cleaning of homo restraints
[unres.git] / source / cluster / wham / src-M / readrtns.F
1       subroutine read_control
2 C
3 C Read molecular data
4 C
5       implicit none
6       include 'DIMENSIONS'
7       include 'sizesclu.dat'
8       include 'COMMON.IOUNITS'
9       include 'COMMON.TIME1'
10       include 'COMMON.SBRIDGE'
11       include 'COMMON.CONTROL'
12       include 'COMMON.CLUSTER'
13       include 'COMMON.CHAIN'
14       include 'COMMON.HEADER'
15       include 'COMMON.FFIELD'
16       include 'COMMON.FREE'
17       include 'COMMON.INTERACT'
18       include "COMMON.SPLITELE"
19       include 'COMMON.SHIELD'
20       character*320 controlcard,ucase
21 #ifdef MPL
22       include 'COMMON.INFO'
23 #endif
24       integer i,i1,i2,it1,it2
25       double precision pi
26       read (INP,'(a80)') titel
27       call card_concat(controlcard)
28
29       call readi(controlcard,'NRES',nres,0)
30       call readi(controlcard,'RESCALE',rescale_mode,2)
31       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
32       write (iout,*) "DISTCHAINMAX",distchainmax
33 C Reading the dimensions of box in x,y,z coordinates
34       call reada(controlcard,'BOXX',boxxsize,100.0d0)
35       call reada(controlcard,'BOXY',boxysize,100.0d0)
36       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
37 c Cutoff range for interactions
38       call reada(controlcard,"R_CUT",r_cut,15.0d0)
39       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
40       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
41       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
42       if (lipthick.gt.0.0d0) then
43        bordliptop=(boxzsize+lipthick)/2.0
44        bordlipbot=bordliptop-lipthick
45 C      endif
46       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
47      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
48       buflipbot=bordlipbot+lipbufthick
49       bufliptop=bordliptop-lipbufthick
50       if ((lipbufthick*2.0d0).gt.lipthick)
51      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
52       endif
53       write(iout,*) "bordliptop=",bordliptop
54       write(iout,*) "bordlipbot=",bordlipbot
55       write(iout,*) "bufliptop=",bufliptop
56       write(iout,*) "buflipbot=",buflipbot
57 C Shielding mode
58       call readi(controlcard,'SHIELD',shield_mode,0)
59       write (iout,*) "SHIELD MODE",shield_mode
60       if (shield_mode.gt.0) then
61       pi=3.141592d0
62 C VSolvSphere the volume of solving sphere
63 C      print *,pi,"pi"
64 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
65 C there will be no distinction between proline peptide group and normal peptide
66 C group in case of shielding parameters
67       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
68       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
69       write (iout,*) VSolvSphere,VSolvSphere_div
70 C long axis of side chain 
71       do i=1,ntyp
72       long_r_sidechain(i)=vbldsc0(1,i)
73       short_r_sidechain(i)=sigma0(i)
74       enddo
75       buff_shield=1.0d0
76       endif
77       call readi(controlcard,'PDBOUT',outpdb,0)
78       call readi(controlcard,'MOL2OUT',outmol2,0)
79       refstr=(index(controlcard,'REFSTR').gt.0)
80       write (iout,*) "REFSTR",refstr
81       pdbref=(index(controlcard,'PDBREF').gt.0)
82       iscode=index(controlcard,'ONE_LETTER')
83       tree=(index(controlcard,'MAKE_TREE').gt.0)
84       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
85       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
86       write (iout,*) "with_dihed_constr ",with_dihed_constr,
87      & " CONSTR_DIST",constr_dist
88       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
89       write (iout,*) "with_theta_constr ",with_theta_constr
90       call flush(iout)
91       min_var=(index(controlcard,'MINVAR').gt.0)
92       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
93       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
94       call readi(controlcard,'NCUT',ncut,1)
95       call readi(controlcard,'SYM',symetr,1)
96       write (iout,*) 'sym', symetr
97       call readi(controlcard,'NSTART',nstart,0)
98       call readi(controlcard,'NEND',nend,0)
99       call reada(controlcard,'ECUT',ecut,10.0d0)
100       call reada(controlcard,'PROB',prob_limit,0.99d0)
101       write (iout,*) "Probability limit",prob_limit
102       lgrp=(index(controlcard,'LGRP').gt.0)
103       caonly=(index(controlcard,'CA_ONLY').gt.0)
104       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
105       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
106       call readi(controlcard,'IOPT',iopt,2) 
107       lside = index(controlcard,"SIDE").gt.0
108       efree = index(controlcard,"EFREE").gt.0
109       call readi(controlcard,'NTEMP',nT,1)
110       write (iout,*) "nT",nT
111       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
112       write (iout,*) "nT",nT
113       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
114       do i=1,nT
115         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
116       enddo
117       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
118       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
119       lprint_int=index(controlcard,"PRINT_INT") .gt.0
120       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
121       write (iout,*) "with_homology_constr ",with_dihed_constr,
122      & " CONSTR_HOMOLOGY",constr_homology
123       print_homology_restraints=
124      & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
125       print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
126       print_homology_models=
127      & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
128       if (min_var) iopt=1
129       return
130       end
131 c--------------------------------------------------------------------------
132       subroutine molread
133 C
134 C Read molecular data.
135 C
136       implicit none
137       include 'DIMENSIONS'
138       include 'COMMON.IOUNITS'
139       include 'COMMON.GEO'
140       include 'COMMON.VAR'
141       include 'COMMON.INTERACT'
142       include 'COMMON.LOCAL'
143       include 'COMMON.NAMES'
144       include 'COMMON.CHAIN'
145       include 'COMMON.FFIELD'
146       include 'COMMON.SBRIDGE'
147       include 'COMMON.HEADER'
148       include 'COMMON.CONTROL'
149       include 'COMMON.CONTACTS'
150       include 'COMMON.TIME1'
151       include 'COMMON.TORCNSTR'
152       include 'COMMON.SHIELD'
153 #ifdef MPL
154       include 'COMMON.INFO'
155 #endif
156       character*4 sequence(maxres)
157       character*800 weightcard
158       integer rescode
159       double precision x(maxvar)
160       integer itype_pdb(maxres)
161       logical seq_comp
162       integer i,j,kkk,i1,i2,it1,it2
163 C
164 C Body
165 C
166 C Read weights of the subsequent energy terms.
167       call card_concat(weightcard)
168       write(iout,*) weightcard
169 C      call reada(weightcard,'WSC',wsc,1.0d0)
170       write(iout,*) wsc
171       call reada(weightcard,'WLONG',wsc,wsc)
172       call reada(weightcard,'WSCP',wscp,1.0d0)
173       call reada(weightcard,'WELEC',welec,1.0D0)
174       call reada(weightcard,'WVDWPP',wvdwpp,welec)
175       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
176       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
177       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
178       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
179       call reada(weightcard,'WTURN3',wturn3,1.0D0)
180       call reada(weightcard,'WTURN4',wturn4,1.0D0)
181       call reada(weightcard,'WTURN6',wturn6,1.0D0)
182       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
183       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
184       call reada(weightcard,'WBOND',wbond,1.0D0)
185       call reada(weightcard,'WTOR',wtor,1.0D0)
186       call reada(weightcard,'WTORD',wtor_d,1.0D0)
187       call reada(weightcard,'WANG',wang,1.0D0)
188       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
189       call reada(weightcard,'SCAL14',scal14,0.4D0)
190       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
191       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
192       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
193       if (index(weightcard,'SOFT').gt.0) ipot=6
194       call reada(weightcard,"D0CM",d0cm,3.78d0)
195       call reada(weightcard,"AKCM",akcm,15.1d0)
196       call reada(weightcard,"AKTH",akth,11.0d0)
197       call reada(weightcard,"AKCT",akct,12.0d0)
198       call reada(weightcard,"V1SS",v1ss,-1.08d0)
199       call reada(weightcard,"V2SS",v2ss,7.61d0)
200       call reada(weightcard,"V3SS",v3ss,13.7d0)
201       call reada(weightcard,"EBR",ebr,-5.50D0)
202       call reada(weightcard,'WSHIELD',wshield,1.0d0)
203       write(iout,*) 'WSHIELD',wshield
204        call reada(weightcard,'WLT',wliptran,0.0D0)
205       call reada(weightcard,"ATRISS",atriss,0.301D0)
206       call reada(weightcard,"BTRISS",btriss,0.021D0)
207       call reada(weightcard,"CTRISS",ctriss,1.001D0)
208       call reada(weightcard,"DTRISS",dtriss,1.001D0)
209       write (iout,*) "ATRISS=", atriss
210       write (iout,*) "BTRISS=", btriss
211       write (iout,*) "CTRISS=", ctriss
212       write (iout,*) "DTRISS=", dtriss
213       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
214       do i=1,maxres
215         dyn_ss_mask(i)=.false.
216       enddo
217       do i=1,maxres-1
218         do j=i+1,maxres
219           dyn_ssbond_ij(i,j)=1.0d300
220         enddo
221       enddo
222       call reada(weightcard,"HT",Ht,0.0D0)
223       if (dyn_ss) then
224         ss_depth=ebr/wsc-0.25*eps(1,1)
225         Ht=Ht/wsc-0.25*eps(1,1)
226         akcm=akcm*wstrain/wsc
227         akth=akth*wstrain/wsc
228         akct=akct*wstrain/wsc
229         v1ss=v1ss*wstrain/wsc
230         v2ss=v2ss*wstrain/wsc
231         v3ss=v3ss*wstrain/wsc
232       else
233         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
234       endif
235       write (iout,'(/a)') "Disulfide bridge parameters:"
236       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
237       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
238       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
239       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
240      & ' v3ss:',v3ss
241
242 C 12/1/95 Added weight for the multi-body term WCORR
243       call reada(weightcard,'WCORRH',wcorr,1.0D0)
244       if (wcorr4.gt.0.0d0) wcorr=wcorr4
245       weights(1)=wsc
246       weights(2)=wscp
247       weights(3)=welec
248       weights(4)=wcorr
249       weights(5)=wcorr5
250       weights(6)=wcorr6
251       weights(7)=wel_loc
252       weights(8)=wturn3
253       weights(9)=wturn4
254       weights(10)=wturn6
255       weights(11)=wang
256       weights(12)=wscloc
257       weights(13)=wtor
258       weights(14)=wtor_d
259       weights(15)=wstrain
260       weights(16)=wvdwpp
261       weights(17)=wbond
262       weights(18)=scal14
263       weights(19)=wsccor
264       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
265      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
266      &  wturn4,wturn6,wsccor
267    10 format (/'Energy-term weights (unscaled):'//
268      & 'WSCC=   ',f10.6,' (SC-SC)'/
269      & 'WSCP=   ',f10.6,' (SC-p)'/
270      & 'WELEC=  ',f10.6,' (p-p electr)'/
271      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
272      & 'WBOND=  ',f10.6,' (stretching)'/
273      & 'WANG=   ',f10.6,' (bending)'/
274      & 'WSCLOC= ',f10.6,' (SC local)'/
275      & 'WTOR=   ',f10.6,' (torsional)'/
276      & 'WTORD=  ',f10.6,' (double torsional)'/
277      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
278      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
279      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
280      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
281      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
282      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
283      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
284      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
285      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
286
287       if (wcorr4.gt.0.0d0) then
288         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
289      &   'between contact pairs of peptide groups'
290         write (iout,'(2(a,f5.3/))')
291      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
292      &  'Range of quenching the correlation terms:',2*delt_corr
293       else if (wcorr.gt.0.0d0) then
294         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
295      &   'between contact pairs of peptide groups'
296       endif
297       write (iout,'(a,f8.3)')
298      &  'Scaling factor of 1,4 SC-p interactions:',scal14
299       write (iout,'(a,f8.3)')
300      &  'General scaling factor of SC-p interactions:',scalscp
301       r0_corr=cutoff_corr-delt_corr
302       do i=1,20
303         aad(i,1)=scalscp*aad(i,1)
304         aad(i,2)=scalscp*aad(i,2)
305         bad(i,1)=scalscp*bad(i,1)
306         bad(i,2)=scalscp*bad(i,2)
307       enddo
308
309       call flush(iout)
310       print *,'indpdb=',indpdb,' pdbref=',pdbref
311
312 C Read sequence if not taken from the pdb file.
313       if (iscode.gt.0) then
314         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
315       else
316         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
317       endif
318 C Convert sequence to numeric code
319       do i=1,nres
320         itype(i)=rescode(i,sequence(i),iscode)
321       enddo
322       print *,nres
323       print '(20i4)',(itype(i),i=1,nres)
324
325       do i=1,nres
326 #ifdef PROCOR
327         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
328 #else
329         if (itype(i).eq.ntyp1) then
330 #endif
331           itel(i)=0
332 #ifdef PROCOR
333         else if (iabs(itype(i+1)).ne.20) then
334 #else
335         else if (iabs(itype(i)).ne.20) then
336 #endif
337           itel(i)=1
338         else
339           itel(i)=2
340         endif
341       enddo
342       write (iout,*) "ITEL"
343       do i=1,nres-1
344         write (iout,*) i,itype(i),itel(i)
345       enddo
346
347       print *,'Call Read_Bridge.'
348       call read_bridge
349 C this fragment reads diheadral constrains
350       if (with_dihed_constr) then
351
352       read (inp,*) ndih_constr
353       if (ndih_constr.gt.0) then
354 C        read (inp,*) ftors
355 C        write (iout,*) 'FTORS',ftors
356 C ftors is the force constant for torsional quartic constrains
357         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
358      &                i=1,ndih_constr)
359         write (iout,*)
360      &   'There are',ndih_constr,' constraints on phi angles.'
361         do i=1,ndih_constr
362           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
363      &  ftors(i)
364         enddo
365         do i=1,ndih_constr
366           phi0(i)=deg2rad*phi0(i)
367           drange(i)=deg2rad*drange(i)
368         enddo
369       endif ! endif ndif_constr.gt.0
370       endif ! with_dihed_constr
371       if (with_theta_constr) then
372 C with_theta_constr is keyword allowing for occurance of theta constrains
373       read (inp,*) ntheta_constr
374 C ntheta_constr is the number of theta constrains
375       if (ntheta_constr.gt.0) then
376 C        read (inp,*) ftors
377         read (inp,*) (itheta_constr(i),theta_constr0(i),
378      &  theta_drange(i),for_thet_constr(i),
379      &  i=1,ntheta_constr)
380 C the above code reads from 1 to ntheta_constr 
381 C itheta_constr(i) residue i for which is theta_constr
382 C theta_constr0 the global minimum value
383 C theta_drange is range for which there is no energy penalty
384 C for_thet_constr is the force constant for quartic energy penalty
385 C E=k*x**4 
386 C        if(me.eq.king.or..not.out1file)then
387          write (iout,*)
388      &   'There are',ntheta_constr,' constraints on phi angles.'
389          do i=1,ntheta_constr
390           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
391      &    theta_drange(i),
392      &    for_thet_constr(i)
393          enddo
394 C        endif
395         do i=1,ntheta_constr
396           theta_constr0(i)=deg2rad*theta_constr0(i)
397           theta_drange(i)=deg2rad*theta_drange(i)
398         enddo
399 C        if(me.eq.king.or..not.out1file)
400 C     &   write (iout,*) 'FTORS',ftors
401 C        do i=1,ntheta_constr
402 C          ii = itheta_constr(i)
403 C          thetabound(1,ii) = phi0(i)-drange(i)
404 C          thetabound(2,ii) = phi0(i)+drange(i)
405 C        enddo
406       endif ! ntheta_constr.gt.0
407       endif! with_theta_constr
408
409       nnt=1
410       nct=nres
411       print *,'NNT=',NNT,' NCT=',NCT
412       if (itype(1).eq.ntyp1) nnt=2
413       if (itype(nres).eq.ntyp1) nct=nct-1
414       if (nstart.lt.nnt) nstart=nnt
415       if (nend.gt.nct .or. nend.eq.0) nend=nct
416       write (iout,*) "nstart",nstart," nend",nend
417       nres0=nres
418       if (constr_homology.gt.0) then
419         call read_constr_homology(print_homology_restraints)
420       endif
421
422 c      if (pdbref) then
423 c        read(inp,'(a)') pdbfile
424 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
425 c        open(ipdbin,file=pdbfile,status='old',err=33)
426 c        goto 34 
427 c  33    write (iout,'(a)') 'Error opening PDB file.'
428 c        stop
429 c  34    continue
430 c        print *,'Begin reading pdb data'
431 c        call readpdb
432 c        print *,'Finished reading pdb data'
433 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
434 c        do i=1,nres
435 c          itype_pdb(i)=itype(i)
436 c        enddo
437 c        close (ipdbin)
438 c        write (iout,'(a,i3)') 'nsup=',nsup
439 c        nstart_seq=nnt
440 c        if (nsup.le.(nct-nnt+1)) then
441 c          do i=0,nct-nnt+1-nsup
442 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
443 c              nstart_seq=nnt+i
444 c              goto 111
445 c            endif
446 c          enddo
447 c          write (iout,'(a)') 
448 c     &            'Error - sequences to be superposed do not match.'
449 c          stop
450 c        else
451 c          do i=0,nsup-(nct-nnt+1)
452 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
453 c     &      then
454 c              nstart_sup=nstart_sup+i
455 c              nsup=nct-nnt+1
456 c              goto 111
457 c            endif
458 c          enddo 
459 c          write (iout,'(a)') 
460 c     &            'Error - sequences to be superposed do not match.'
461 c        endif
462 c  111   continue
463 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
464 c     &                 ' nstart_seq=',nstart_seq
465 c      endif
466       call init_int_table
467       call setup_var
468       write (iout,*) "molread: REFSTR",refstr
469       if (refstr) then
470         if (.not.pdbref) then
471           call read_angles(inp,*38)
472           goto 39
473    38     write (iout,'(a)') 'Error reading reference structure.'
474 #ifdef MPL
475           call mp_stopall(Error_Msg)
476 #else
477           stop 'Error reading reference structure'
478 #endif
479    39     call chainbuild     
480           nstart_sup=nnt
481           nstart_seq=nnt
482           nsup=nct-nnt+1
483           kkk=1
484           do i=1,2*nres
485             do j=1,3
486               cref(j,i,kkk)=c(j,i)
487             enddo
488           enddo
489         endif
490         call contact(.true.,ncont_ref,icont_ref)
491       endif
492        if (ns.gt.0) then
493 C        write (iout,'(/a,i3,a)')
494 C     &  'The chain contains',ns,' disulfide-bridging cysteines.'
495         write (iout,'(20i4)') (iss(i),i=1,ns)
496        if (dyn_ss) then
497           write(iout,*)"Running with dynamic disulfide-bond formation"
498        else
499         write (iout,'(/a/)') 'Pre-formed links are:'
500         do i=1,nss
501           i1=ihpb(i)-nres
502           i2=jhpb(i)-nres
503           it1=itype(i1)
504           it2=itype(i2)
505           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
506      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
507      &    ebr,forcon(i)
508         enddo
509         write (iout,'(a)')
510        endif
511       endif
512       if (ns.gt.0.and.dyn_ss) then
513           do i=nss+1,nhpb
514             ihpb(i-nss)=ihpb(i)
515             jhpb(i-nss)=jhpb(i)
516             forcon(i-nss)=forcon(i)
517             dhpb(i-nss)=dhpb(i)
518           enddo
519           nhpb=nhpb-nss
520           nss=0
521           call hpb_partition
522           do i=1,ns
523             dyn_ss_mask(iss(i))=.true.
524           enddo
525       endif
526 c Read distance restraints
527       if (constr_dist.gt.0) then
528         call read_dist_constr
529         call hpb_partition
530       endif
531       return
532       end
533 c-----------------------------------------------------------------------------
534       logical function seq_comp(itypea,itypeb,length)
535       implicit none
536       integer length,itypea(length),itypeb(length)
537       integer i
538       do i=1,length
539         if (itypea(i).ne.itypeb(i)) then
540           seq_comp=.false.
541           return
542         endif
543       enddo
544       seq_comp=.true.
545       return
546       end
547 c-----------------------------------------------------------------------------
548       subroutine read_bridge
549 C Read information about disulfide bridges.
550       implicit none
551       include 'DIMENSIONS'
552       include 'COMMON.IOUNITS'
553       include 'COMMON.GEO'
554       include 'COMMON.VAR'
555       include 'COMMON.INTERACT'
556       include 'COMMON.LOCAL'
557       include 'COMMON.NAMES'
558       include 'COMMON.CHAIN'
559       include 'COMMON.FFIELD'
560       include 'COMMON.SBRIDGE'
561       include 'COMMON.HEADER'
562       include 'COMMON.CONTROL'
563       include 'COMMON.TIME1'
564 #ifdef MPL
565       include 'COMMON.INFO'
566 #endif
567       integer i,j
568 C Read bridging residues.
569       read (inp,*) ns,(iss(i),i=1,ns)
570       print *,'ns=',ns
571 C Check whether the specified bridging residues are cystines.
572       do i=1,ns
573         if (itype(iss(i)).ne.1) then
574           write (iout,'(2a,i3,a)') 
575      &   'Do you REALLY think that the residue ',
576      &    restyp(itype(iss(i))),i,
577      &   ' can form a disulfide bridge?!!!'
578           write (*,'(2a,i3,a)') 
579      &   'Do you REALLY think that the residue ',
580      &   restyp(itype(iss(i))),i,
581      &   ' can form a disulfide bridge?!!!'
582 #ifdef MPL
583          call mp_stopall(error_msg)
584 #else
585          stop
586 #endif
587         endif
588       enddo
589 C Read preformed bridges.
590       if (ns.gt.0) then
591       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
592       if (nss.gt.0) then
593         nhpb=nss
594 C Check if the residues involved in bridges are in the specified list of
595 C bridging residues.
596         do i=1,nss
597           do j=1,i-1
598             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
599      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
600               write (iout,'(a,i3,a)') 'Disulfide pair',i,
601      &      ' contains residues present in other pairs.'
602               write (*,'(a,i3,a)') 'Disulfide pair',i,
603      &      ' contains residues present in other pairs.'
604 #ifdef MPL
605               call mp_stopall(error_msg)
606 #else
607               stop 
608 #endif
609             endif
610           enddo
611           do j=1,ns
612             if (ihpb(i).eq.iss(j)) goto 10
613           enddo
614           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
615    10     continue
616           do j=1,ns
617             if (jhpb(i).eq.iss(j)) goto 20
618           enddo
619           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
620    20     continue
621 C          dhpb(i)=dbr
622 C          forcon(i)=fbr
623         enddo
624         do i=1,nss
625           ihpb(i)=ihpb(i)+nres
626           jhpb(i)=jhpb(i)+nres
627         enddo
628       endif
629       endif
630       return
631       end
632 c----------------------------------------------------------------------------
633       subroutine read_angles(kanal,*)
634       implicit none
635       include 'DIMENSIONS'
636       include 'COMMON.GEO'
637       include 'COMMON.VAR'
638       include 'COMMON.CHAIN'
639       include 'COMMON.IOUNITS'
640       integer i,kanal
641       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
642       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
643       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
644       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
645       do i=1,nres
646         theta(i)=deg2rad*theta(i)
647         phi(i)=deg2rad*phi(i)
648         alph(i)=deg2rad*alph(i)
649         omeg(i)=deg2rad*omeg(i)
650       enddo
651       return
652    10 return1
653       end
654 c----------------------------------------------------------------------------
655       subroutine reada(rekord,lancuch,wartosc,default)
656       implicit none
657       character*(*) rekord,lancuch
658       double precision wartosc,default
659       integer ilen,iread
660       external ilen
661       iread=index(rekord,lancuch)
662       if (iread.eq.0) then
663         wartosc=default 
664         return
665       endif   
666       iread=iread+ilen(lancuch)+1
667       read (rekord(iread:),*) wartosc
668       return
669       end
670 c----------------------------------------------------------------------------
671       subroutine multreada(rekord,lancuch,tablica,dim,default)
672       implicit none
673       integer dim,i
674       double precision tablica(dim),default
675       character*(*) rekord,lancuch
676       integer ilen,iread
677       external ilen
678       do i=1,dim
679         tablica(i)=default 
680       enddo
681       iread=index(rekord,lancuch)
682       if (iread.eq.0) return
683       iread=iread+ilen(lancuch)+1
684       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
685    10 return
686       end
687 c----------------------------------------------------------------------------
688       subroutine readi(rekord,lancuch,wartosc,default)
689       implicit none
690       character*(*) rekord,lancuch
691       integer wartosc,default
692       integer ilen,iread
693       external ilen
694       iread=index(rekord,lancuch)
695       if (iread.eq.0) then
696         wartosc=default 
697         return
698       endif   
699       iread=iread+ilen(lancuch)+1
700       read (rekord(iread:),*) wartosc
701       return
702       end
703 C----------------------------------------------------------------------
704       subroutine multreadi(rekord,lancuch,tablica,dim,default)
705       implicit none
706       integer dim,i
707       integer tablica(dim),default
708       character*(*) rekord,lancuch
709       character*80 aux
710       integer ilen,iread
711       external ilen
712       do i=1,dim
713         tablica(i)=default
714       enddo
715       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
716       if (iread.eq.0) return
717       iread=iread+ilen(lancuch)+1
718       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
719    10 return
720       end
721
722 c----------------------------------------------------------------------------
723       subroutine card_concat(card)
724       include 'DIMENSIONS'
725       include 'COMMON.IOUNITS'
726       character*(*) card
727       character*80 karta,ucase
728       external ilen
729       read (inp,'(a)') karta
730       karta=ucase(karta)
731       card=' '
732       do while (karta(80:80).eq.'&')
733         card=card(:ilen(card)+1)//karta(:79)
734         read (inp,'(a)') karta
735         karta=ucase(karta)
736       enddo
737       card=card(:ilen(card)+1)//karta
738       return
739       end
740 c----------------------------------------------------------------------------
741       subroutine openunits
742       implicit none
743       include 'DIMENSIONS'    
744 #ifdef MPI
745       include "mpif.h"
746       character*3 liczba
747       include "COMMON.MPI"
748 #endif
749       include 'COMMON.IOUNITS'
750       include 'COMMON.CONTROL'
751       integer lenpre,lenpot,ilen
752       external ilen
753       character*16 cformat,cprint
754       character*16 ucase
755       integer lenint,lenout
756       call getenv('INPUT',prefix)
757       call getenv('OUTPUT',prefout)
758       call getenv('INTIN',prefintin)
759       call getenv('COORD',cformat)
760       call getenv('PRINTCOOR',cprint)
761       call getenv('SCRATCHDIR',scratchdir)
762       from_bx=.true.
763       from_cx=.false.
764       if (index(ucase(cformat),'CX').gt.0) then
765         from_cx=.true.
766         from_bx=.false.
767       endif
768       from_cart=.true.
769       lenpre=ilen(prefix)
770       lenout=ilen(prefout)
771       lenint=ilen(prefintin)
772 C Get the names and open the input files
773       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
774 #ifdef MPI
775       write (liczba,'(bz,i3.3)') me
776       outname=prefout(:lenout)//'_clust.out_'//liczba
777 #else
778       outname=prefout(:lenout)//'_clust.out'
779 #endif
780       if (from_bx) then
781         intinname=prefintin(:lenint)//'.bx'
782       else if (from_cx) then
783         intinname=prefintin(:lenint)//'.cx'
784       else
785         intinname=prefintin(:lenint)//'.int'
786       endif
787       rmsname=prefintin(:lenint)//'.rms'
788       open (jplot,file=prefout(:ilen(prefout))//'.tex',
789      &       status='unknown')
790       open (jrms,file=rmsname,status='unknown')
791       open(iout,file=outname,status='unknown')
792 C Get parameter filenames and open the parameter files.
793       call getenv('BONDPAR',bondname)
794       open (ibond,file=bondname,status='old')
795       call getenv('THETPAR',thetname)
796       open (ithep,file=thetname,status='old')
797       call getenv('ROTPAR',rotname)
798       open (irotam,file=rotname,status='old')
799       call getenv('TORPAR',torname)
800       open (itorp,file=torname,status='old')
801       call getenv('TORDPAR',tordname)
802       open (itordp,file=tordname,status='old')
803       call getenv('FOURIER',fouriername)
804       open (ifourier,file=fouriername,status='old')
805       call getenv('ELEPAR',elename)
806       open (ielep,file=elename,status='old')
807       call getenv('SIDEPAR',sidename)
808       open (isidep,file=sidename,status='old')
809       call getenv('SIDEP',sidepname)
810       open (isidep1,file=sidepname,status="old")
811       call getenv('SCCORPAR',sccorname)
812       open (isccor,file=sccorname,status="old")
813       call getenv('LIPTRANPAR',liptranname)
814       open (iliptranpar,file=liptranname,status='old')
815 #ifndef OLDSCP
816 C
817 C 8/9/01 In the newest version SCp interaction constants are read from a file
818 C Use -DOLDSCP to use hard-coded constants instead.
819 C
820       call getenv('SCPPAR',scpname)
821       open (iscpp,file=scpname,status='old')
822 #endif
823       return
824       end
825       subroutine read_dist_constr
826       implicit real*8 (a-h,o-z)
827       include 'DIMENSIONS'
828 #ifdef MPI
829       include 'mpif.h'
830 #endif
831       include 'COMMON.CONTROL'
832       include 'COMMON.CHAIN'
833       include 'COMMON.IOUNITS'
834       include 'COMMON.SBRIDGE'
835       integer ifrag_(2,100),ipair_(2,100)
836       double precision wfrag_(100),wpair_(100)
837       character*500 controlcard
838       logical lprn /.true./
839       write (iout,*) "Calling read_dist_constr"
840 C      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
841 C      call flush(iout)
842       write(iout,*) "TU sie wywalam?"
843       call card_concat(controlcard)
844       write (iout,*) controlcard
845       call flush(iout)
846       call readi(controlcard,"NFRAG",nfrag_,0)
847       call readi(controlcard,"NPAIR",npair_,0)
848       call readi(controlcard,"NDIST",ndist_,0)
849       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
850       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
851       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
852       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
853       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
854       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
855       write (iout,*) "IFRAG"
856       do i=1,nfrag_
857         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
858       enddo
859       write (iout,*) "IPAIR"
860       do i=1,npair_
861         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
862       enddo
863       call flush(iout)
864       do i=1,nfrag_
865         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
866         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
867      &    ifrag_(2,i)=nstart_sup+nsup-1
868 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
869         call flush(iout)
870         if (wfrag_(i).gt.0.0d0) then
871         do j=ifrag_(1,i),ifrag_(2,i)-1
872           do k=j+1,ifrag_(2,i)
873             write (iout,*) "j",j," k",k
874             ddjk=dist(j,k)
875             if (constr_dist.eq.1) then
876               nhpb=nhpb+1
877               ihpb(nhpb)=j
878               jhpb(nhpb)=k
879               dhpb(nhpb)=ddjk
880               forcon(nhpb)=wfrag_(i) 
881             else if (constr_dist.eq.2) then
882               if (ddjk.le.dist_cut) then
883                 nhpb=nhpb+1
884                 ihpb(nhpb)=j
885                 jhpb(nhpb)=k
886                 dhpb(nhpb)=ddjk
887                 forcon(nhpb)=wfrag_(i) 
888               endif
889             else
890               nhpb=nhpb+1
891               ihpb(nhpb)=j
892               jhpb(nhpb)=k
893               dhpb(nhpb)=ddjk
894               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
895             endif
896             if (lprn)
897      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
898      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
899           enddo
900         enddo
901         endif
902       enddo
903       do i=1,npair_
904         if (wpair_(i).gt.0.0d0) then
905         ii = ipair_(1,i)
906         jj = ipair_(2,i)
907         if (ii.gt.jj) then
908           itemp=ii
909           ii=jj
910           jj=itemp
911         endif
912         do j=ifrag_(1,ii),ifrag_(2,ii)
913           do k=ifrag_(1,jj),ifrag_(2,jj)
914             nhpb=nhpb+1
915             ihpb(nhpb)=j
916             jhpb(nhpb)=k
917             forcon(nhpb)=wpair_(i)
918             dhpb(nhpb)=dist(j,k)
919             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
920      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
921           enddo
922         enddo
923         endif
924       enddo 
925       do i=1,ndist_
926         if (constr_dist.eq.11) then
927         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
928      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
929         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
930 C        write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
931 C     &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
932         else
933         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
934         endif
935         if (forcon(nhpb+1).gt.0.0d0) then
936           nhpb=nhpb+1
937           if (ibecarb(i).gt.0) then
938             ihpb(i)=ihpb(i)+nres
939             jhpb(i)=jhpb(i)+nres
940           endif
941           if (dhpb(nhpb).eq.0.0d0)
942      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
943 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
944           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
945      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
946         endif
947 C      endif
948       enddo
949       call hpb_partition
950 #ifdef AIX
951       call flush_(iout)
952 #else
953       call flush(iout)
954 #endif
955       return
956       end
957
958 c====-------------------------------------------------------------------
959       subroutine read_constr_homology(lprn)
960
961       include 'DIMENSIONS'
962 #ifdef MPI
963       include 'mpif.h'
964 #endif
965       include 'COMMON.SETUP'
966       include 'COMMON.CONTROL'
967       include 'COMMON.CHAIN'
968       include 'COMMON.IOUNITS'
969       include 'COMMON.GEO'
970       include 'COMMON.INTERACT'
971       include 'COMMON.NAMES'
972       include 'COMMON.HOMRESTR'
973 c
974 c For new homol impl
975 c
976       include 'COMMON.VAR'
977 c     include 'include_unres/COMMON.VAR'
978 c
979
980 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
981 c    &                 dist_cut
982 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
983 c    &    sigma_odl_temp(maxres,maxres,max_template)
984       character*2 kic2
985       character*24 model_ki_dist, model_ki_angle
986       character*500 controlcard
987       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
988       integer idomain(max_template,maxres)
989       integer ilen
990       external ilen
991       logical lprn
992       logical unres_pdb,liiflag
993 c
994 c     FP - Nov. 2014 Temporary specifications for new vars
995 c
996       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
997      &                 rescore3_tmp
998       double precision, dimension (max_template,maxres) :: rescore
999       double precision, dimension (max_template,maxres) :: rescore2
1000       double precision, dimension (max_template,maxres) :: rescore3
1001       character*24 tpl_k_rescore
1002 c -----------------------------------------------------------------
1003 c Reading multiple PDB ref structures and calculation of retraints
1004 c not using pre-computed ones stored in files model_ki_{dist,angle}
1005 c FP (Nov., 2014)
1006 c -----------------------------------------------------------------
1007 c
1008 c
1009 c Alternative: reading from input
1010 #ifdef DEBUG
1011       write (iout,*) "BEGIN READ HOMOLOGY INFO"
1012 #ifdef AIX
1013       call flush_(iout)
1014 #else
1015       call flush(iout)
1016 #endif
1017 #endif
1018       call card_concat(controlcard)
1019       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1020       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1021       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1022       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1023       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1024       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1025       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1026       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1027       if (homol_nset.gt.1)then
1028          call readi(controlcard,"ISET",iset,1)
1029          call card_concat(controlcard)
1030          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1031       else
1032         iset=1
1033         waga_homology(1)=1.0
1034       endif
1035 c
1036 #ifdef DEBUG
1037       write(iout,*) "read_constr_homology iset",iset
1038       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1039 #ifdef AIX
1040       call flush_(iout)
1041 #else
1042       call flush(iout)
1043 #endif
1044 #endif
1045 cd      write (iout,*) "nnt",nnt," nct",nct
1046 cd      call flush(iout)
1047
1048
1049       lim_odl=0
1050       lim_dih=0
1051 c
1052 c  New
1053 c
1054 c
1055 c  Reading HM global scores (prob not required)
1056 c
1057       do i = nnt,nct
1058         do k=1,constr_homology
1059          idomain(k,i)=0
1060         enddo
1061       enddo
1062 c     open (4,file="HMscore")
1063 c     do k=1,constr_homology
1064 c       read (4,*,end=521) hmscore_tmp
1065 c       hmscore(k)=hmscore_tmp ! Another transformation can be used 
1066 c       write(*,*) "Model", k, ":", hmscore(k)
1067 c     enddo
1068 c521  continue
1069
1070       ii=0
1071       do i = nnt,nct-2 
1072         do j=i+2,nct 
1073         ii=ii+1
1074         ii_in_use(ii)=0
1075         enddo
1076       enddo
1077 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1078
1079       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1080       do k=1,constr_homology
1081
1082         read(inp,'(a)') pdbfile
1083 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1084 c  Next stament causes error upon compilation (?)
1085 c       if(me.eq.king.or. .not. out1file)
1086 c         write (iout,'(2a)') 'PDB data will be read from file ',
1087 c    &   pdbfile(:ilen(pdbfile))
1088          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1089      &  pdbfile(:ilen(pdbfile))
1090         open(ipdbin,file=pdbfile,status='old',err=33)
1091         goto 34
1092   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1093      &  pdbfile(:ilen(pdbfile))
1094         stop
1095   34    continue
1096 c        print *,'Begin reading pdb data'
1097 c
1098 c Files containing res sim or local scores (former containing sigmas)
1099 c
1100
1101         write(kic2,'(bz,i2.2)') k
1102
1103         tpl_k_rescore="template"//kic2//".sco"
1104
1105         unres_pdb=.false.
1106         call readpdb
1107         do i=1,2*nres
1108           do j=1,3
1109             crefjlee(j,i)=c(j,i)
1110           enddo
1111         enddo
1112 #ifdef DEBUG
1113         do i=1,nres
1114           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1115      &      (crefjlee(j,i+nres),j=1,3)
1116         enddo
1117         write (iout,*) "READ HOMOLOGY INFO"
1118         write (iout,*) "read_constr_homology x: after reading pdb file"
1119         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1120         write (iout,*) "waga_dist",waga_dist
1121         write (iout,*) "waga_angle",waga_angle
1122         write (iout,*) "waga_theta",waga_theta
1123         write (iout,*) "waga_d",waga_d
1124         write (iout,*) "dist_cut",dist_cut
1125 #endif
1126 #ifdef AIX
1127       call flush_(iout)
1128 #else
1129       call flush(iout)
1130 #endif
1131
1132 c
1133 c     Distance restraints
1134 c
1135 c          ... --> odl(k,ii)
1136 C Copy the coordinates from reference coordinates (?)
1137         do i=1,2*nres
1138           do j=1,3
1139 c            c(j,i)=cref(j,i)
1140 c           write (iout,*) "c(",j,i,") =",c(j,i)
1141           enddo
1142         enddo
1143 c
1144 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1145 c
1146 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1147           open (ientin,file=tpl_k_rescore,status='old')
1148           if (nnt.gt.1) rescore(k,1)=0.0d0
1149           do irec=nnt,nct ! loop for reading res sim 
1150             if (read2sigma) then
1151              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1152      &                                rescore3_tmp,idomain_tmp
1153              i_tmp=i_tmp+nnt-1
1154              idomain(k,i_tmp)=idomain_tmp
1155              rescore(k,i_tmp)=rescore_tmp
1156              rescore2(k,i_tmp)=rescore2_tmp
1157              rescore3(k,i_tmp)=rescore3_tmp
1158              write(iout,'(a7,i5,3f10.5,i5)') "rescore",
1159      &                      i_tmp,rescore2_tmp,rescore_tmp,
1160      &                                rescore3_tmp,idomain_tmp
1161             else
1162              idomain(k,irec)=1
1163              read (ientin,*,end=1401) rescore_tmp
1164
1165 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1166              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1167 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1168             endif
1169           enddo  
1170  1401   continue
1171         close (ientin)        
1172         if (waga_dist.ne.0.0d0) then
1173           ii=0
1174           do i = nnt,nct-2 
1175             do j=i+2,nct 
1176
1177               x12=c(1,i)-c(1,j)
1178               y12=c(2,i)-c(2,j)
1179               z12=c(3,i)-c(3,j)
1180               distal=dsqrt(x12*x12+y12*y12+z12*z12) 
1181 c              write (iout,*) k,i,j,distal,dist2_cut
1182
1183             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1184      &            .and. distal.le.dist2_cut ) then
1185
1186               ii=ii+1
1187               ii_in_use(ii)=1
1188               l_homo(k,ii)=.true.
1189
1190 c             write (iout,*) "k",k
1191 c             write (iout,*) "i",i," j",j," constr_homology",
1192 c    &                       constr_homology
1193               ires_homo(ii)=i
1194               jres_homo(ii)=j
1195               odl(k,ii)=distal
1196               if (read2sigma) then
1197                 sigma_odl(k,ii)=0
1198                 do ik=i,j
1199                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1200                 enddo
1201                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1202                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = 
1203      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1204               else
1205                 if (odl(k,ii).le.dist_cut) then
1206                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) 
1207                 else
1208 #ifdef OLDSIGMA
1209                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
1210      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1211 #else
1212                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
1213      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1214 #endif
1215                 endif
1216               endif
1217               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) 
1218             else
1219               ii=ii+1
1220               l_homo(k,ii)=.false.
1221             endif
1222             enddo
1223           enddo
1224         lim_odl=ii
1225         endif
1226 c
1227 c     Theta, dihedral and SC retraints
1228 c
1229         if (waga_angle.gt.0.0d0) then
1230 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1231 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1232 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1233 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1234 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1235 c    &                            sigma_dih(k,i+nnt-1)
1236 c         enddo
1237 c1402   continue
1238 c         close (ientin)
1239           do i = nnt+3,nct 
1240             if (idomain(k,i).eq.0) then 
1241                sigma_dih(k,i)=0.0
1242                cycle
1243             endif
1244             dih(k,i)=phiref(i) ! right?
1245 c           read (ientin,*) sigma_dih(k,i) ! original variant
1246 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1247 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1248 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1249 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1250 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1251
1252             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1253      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1254 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1255 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1256 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1257 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1258 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1259 c   Instead of res sim other local measure of b/b str reliability possible
1260             if (sigma_dih(k,i).ne.0)
1261      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1262 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1263           enddo
1264           lim_dih=nct-nnt-2 
1265         endif
1266
1267         if (waga_theta.gt.0.0d0) then
1268 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1269 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1270 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1271 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1272 c    &                              sigma_theta(k,i+nnt-1)
1273 c         enddo
1274 c1403   continue
1275 c         close (ientin)
1276
1277           do i = nnt+2,nct ! right? without parallel.
1278 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1279 c         do i=ithet_start,ithet_end ! with FG parallel.
1280              if (idomain(k,i).eq.0) then  
1281               sigma_theta(k,i)=0.0
1282               cycle
1283              endif
1284              thetatpl(k,i)=thetaref(i)
1285 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1286 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1287 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1288 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1289 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1290              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1291      &                        rescore(k,i-2))/3.0
1292 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1293              if (sigma_theta(k,i).ne.0)
1294      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1295
1296 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1297 c                             rescore(k,i-2) !  right expression ?
1298 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1299           enddo
1300         endif
1301
1302         if (waga_d.gt.0.0d0) then
1303 c       open (ientin,file=tpl_k_sigma_d,status='old')
1304 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1305 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1306 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1307 c    &                          sigma_d(k,i+nnt-1)
1308 c         enddo
1309 c1404   continue
1310
1311           do i = nnt,nct ! right? without parallel.
1312 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1313 c         do i=loc_start,loc_end ! with FG parallel.
1314                if (itype(i).eq.10) cycle 
1315                if (idomain(k,i).eq.0 ) then 
1316                   sigma_d(k,i)=0.0
1317                   cycle
1318                endif
1319                xxtpl(k,i)=xxref(i)
1320                yytpl(k,i)=yyref(i)
1321                zztpl(k,i)=zzref(i)
1322 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1323 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1324 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1325 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1326                sigma_d(k,i)=rescore3(k,i) !  right expression ?
1327                if (sigma_d(k,i).ne.0)
1328      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1329
1330 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1331 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1332 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1333           enddo
1334         endif
1335       enddo
1336 c
1337 c remove distance restraints not used in any model from the list
1338 c shift data in all arrays
1339 c
1340       if (waga_dist.ne.0.0d0) then
1341         ii=0
1342         liiflag=.true.
1343         do i=nnt,nct-2 
1344          do j=i+2,nct 
1345           ii=ii+1
1346           if (ii_in_use(ii).eq.0.and.liiflag) then
1347             liiflag=.false.
1348             iistart=ii
1349           endif
1350           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1351      &                   .not.liiflag.and.ii.eq.lim_odl) then
1352              if (ii.eq.lim_odl) then
1353               iishift=ii-iistart+1
1354              else
1355               iishift=ii-iistart
1356              endif
1357              liiflag=.true.
1358              do ki=iistart,lim_odl-iishift
1359               ires_homo(ki)=ires_homo(ki+iishift)
1360               jres_homo(ki)=jres_homo(ki+iishift)
1361               ii_in_use(ki)=ii_in_use(ki+iishift)
1362               do k=1,constr_homology
1363                odl(k,ki)=odl(k,ki+iishift)
1364                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1365                l_homo(k,ki)=l_homo(k,ki+iishift)
1366               enddo
1367              enddo
1368              ii=ii-iishift
1369              lim_odl=lim_odl-iishift
1370           endif
1371          enddo
1372         enddo
1373       endif
1374       if (constr_homology.gt.0) call homology_partition
1375       if (constr_homology.gt.0) call init_int_table
1376 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1377 cd     & "lim_xx=",lim_xx
1378 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1379 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1380 c
1381 c Print restraints
1382 c
1383       if (.not.lprn) return
1384 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1385       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1386        write (iout,*) "Distance restraints from templates"
1387        do ii=1,lim_odl
1388        write(iout,'(3i5,100(2f8.2,1x,l1,4x))') 
1389      &  ii,ires_homo(ii),jres_homo(ii),
1390      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1391      &  ki=1,constr_homology)
1392        enddo
1393        write (iout,*) "Dihedral angle restraints from templates"
1394        do i=nnt+3,nct
1395         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1396      &      (rad2deg*dih(ki,i),
1397      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1398        enddo
1399        write (iout,*) "Virtual-bond angle restraints from templates"
1400        do i=nnt+2,nct
1401         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1402      &      (rad2deg*thetatpl(ki,i),
1403      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1404        enddo
1405        write (iout,*) "SC restraints from templates"
1406        do i=nnt,nct
1407         write(iout,'(i5,100(4f8.2,4x))') i,
1408      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1409      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1410        enddo
1411       endif
1412 c -----------------------------------------------------------------
1413       return
1414       end