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