cluster_wham Adam's new constr_dist single chain
[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       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
80
81 #ifdef AIX
82       call flush_(iout)
83 #else
84       call flush(iout)
85 #endif
86       if (min_var) iopt=1
87       return
88       end
89 c--------------------------------------------------------------------------
90       subroutine molread
91 C
92 C Read molecular data.
93 C
94       implicit none
95       include 'DIMENSIONS'
96       include 'COMMON.IOUNITS'
97       include 'COMMON.GEO'
98       include 'COMMON.VAR'
99       include 'COMMON.INTERACT'
100       include 'COMMON.LOCAL'
101       include 'COMMON.NAMES'
102       include 'COMMON.CHAIN'
103       include 'COMMON.FFIELD'
104       include 'COMMON.SBRIDGE'
105       include 'COMMON.HEADER'
106       include 'COMMON.CONTROL'
107       include 'COMMON.CONTACTS'
108       include 'COMMON.TIME1'
109       include 'COMMON.TORCNSTR'
110 #ifdef MPL
111       include 'COMMON.INFO'
112 #endif
113       character*4 sequence(maxres)
114       character*800 weightcard
115       integer rescode
116       double precision x(maxvar)
117       integer itype_pdb(maxres)
118       logical seq_comp
119       integer i,j
120       write (iout,*) " MOLREAD: NRES",NRES
121 C
122 C Body
123 C
124 C Read weights of the subsequent energy terms.
125       call card_concat(weightcard)
126       call reada(weightcard,'WSC',wsc,1.0d0)
127       call reada(weightcard,'WLONG',wsc,wsc)
128       call reada(weightcard,'WSCP',wscp,1.0d0)
129       call reada(weightcard,'WELEC',welec,1.0D0)
130       call reada(weightcard,'WVDWPP',wvdwpp,welec)
131       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
132       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
133       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
134       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
135       call reada(weightcard,'WTURN3',wturn3,1.0D0)
136       call reada(weightcard,'WTURN4',wturn4,1.0D0)
137       call reada(weightcard,'WTURN6',wturn6,1.0D0)
138       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
139       call reada(weightcard,'WBOND',wbond,1.0D0)
140       call reada(weightcard,'WTOR',wtor,1.0D0)
141       call reada(weightcard,'WTORD',wtor_d,1.0D0)
142       call reada(weightcard,'WANG',wang,1.0D0)
143       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
144       call reada(weightcard,'SCAL14',scal14,0.4D0)
145       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
146       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
147       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
148       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
149       call reada(weightcard,"D0CM",d0cm,3.78d0)
150       call reada(weightcard,"AKCM",akcm,15.1d0)
151       call reada(weightcard,"AKTH",akth,11.0d0)
152       call reada(weightcard,"AKCT",akct,12.0d0)
153       call reada(weightcard,"V1SS",v1ss,-1.08d0)
154       call reada(weightcard,"V2SS",v2ss,7.61d0)
155       call reada(weightcard,"V3SS",v3ss,13.7d0)
156       call reada(weightcard,"EBR",ebr,-5.50D0)
157 C     Bartek
158       call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
159       call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
160       call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
161       call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
162       if (index(weightcard,'SOFT').gt.0) ipot=6
163 C 12/1/95 Added weight for the multi-body term WCORR
164       call reada(weightcard,'WCORRH',wcorr,1.0D0)
165       do i=1,maxres-1
166         do j=i+1,maxres
167           dyn_ssbond_ij(i,j)=1.0d300
168         enddo
169       enddo
170       call reada(weightcard,"HT",Ht,0.0D0)
171       if (dyn_ss) then
172         ss_depth=ebr/wsc-0.25*eps(1,1)
173         Ht=Ht/wsc-0.25*eps(1,1)
174         akcm=akcm*wstrain/wsc
175         akth=akth*wstrain/wsc
176         akct=akct*wstrain/wsc
177         v1ss=v1ss*wstrain/wsc
178         v2ss=v2ss*wstrain/wsc
179         v3ss=v3ss*wstrain/wsc
180       else
181         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
182       endif
183       write (iout,'(/a)') "Disulfide bridge parameters:"
184       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
185         write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
186       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
187       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
188       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
189      & ' v3ss:',v3ss
190       write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
191       if (wcorr4.gt.0.0d0) wcorr=wcorr4
192       weights(1)=wsc
193       weights(2)=wscp
194       weights(3)=welec
195       weights(4)=wcorr
196       weights(5)=wcorr5
197       weights(6)=wcorr6
198       weights(7)=wel_loc
199       weights(8)=wturn3
200       weights(9)=wturn4
201       weights(10)=wturn6
202       weights(11)=wang
203       weights(12)=wscloc
204       weights(13)=wtor
205       weights(14)=wtor_d
206       weights(15)=wstrain
207       weights(16)=wvdwpp
208       weights(17)=wbond
209       weights(18)=scal14
210       weights(22)=wdfa_dist
211       weights(23)=wdfa_tor
212       weights(24)=wdfa_nei
213       weights(25)=wdfa_beta
214       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
215      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
216      &  wturn4,wturn6,wsccor,wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
217    10 format (/'Energy-term weights (unscaled):'//
218      & 'WSCC=   ',f10.6,' (SC-SC)'/
219      & 'WSCP=   ',f10.6,' (SC-p)'/
220      & 'WELEC=  ',f10.6,' (p-p electr)'/
221      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
222      & 'WBOND=  ',f10.6,' (stretching)'/
223      & 'WANG=   ',f10.6,' (bending)'/
224      & 'WSCLOC= ',f10.6,' (SC local)'/
225      & 'WTOR=   ',f10.6,' (torsional)'/
226      & 'WTORD=  ',f10.6,' (double torsional)'/
227      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
228      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
229      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
230      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
231      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
232      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
233      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
234      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
235      & 'WSCCOR= ',f10.6,' (SC-backbone torsional correalations)'/
236      & 'WDFAD=  ',f10.6,' (DFA distance)'/
237      & 'WDFAT=  ',f10.6,' (DFA torsional)'/
238      & 'WDFAN=  ',f10.6,' (DFA neighbors)'/
239      & 'WDFAB=  ',f10.6,' (DFA beta)'/)
240       if (wcorr4.gt.0.0d0) then
241         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
242      &   'between contact pairs of peptide groups'
243         write (iout,'(2(a,f5.3/))')
244      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
245      &  'Range of quenching the correlation terms:',2*delt_corr
246       else if (wcorr.gt.0.0d0) then
247         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
248      &   'between contact pairs of peptide groups'
249       endif
250       write (iout,'(a,f8.3)')
251      &  'Scaling factor of 1,4 SC-p interactions:',scal14
252       write (iout,'(a,f8.3)')
253      &  'General scaling factor of SC-p interactions:',scalscp
254       r0_corr=cutoff_corr-delt_corr
255       do i=1,20
256         aad(i,1)=scalscp*aad(i,1)
257         aad(i,2)=scalscp*aad(i,2)
258         bad(i,1)=scalscp*bad(i,1)
259         bad(i,2)=scalscp*bad(i,2)
260       enddo
261
262 #ifdef AIX
263       call flush_(iout)
264 #else
265       call flush(iout)
266 #endif
267       print *,'indpdb=',indpdb,' pdbref=',pdbref
268
269 C Read sequence if not taken from the pdb file.
270       if (iscode.gt.0) then
271         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
272       else
273         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
274       endif
275 C Convert sequence to numeric code
276       do i=1,nres
277         itype(i)=rescode(i,sequence(i),iscode)
278       enddo
279       if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
280         write (iout,*)
281      &   "Glycine is the first full residue, initial dummy deleted"
282         do i=1,nres
283           itype(i)=itype(i+1)
284         enddo
285         nres=nres-1
286       endif
287       if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
288         write (iout,*)
289      &   "Glycine is the last full residue, terminal dummy deleted"
290         nres=nres-1
291       endif
292       print *,nres
293       print '(20i4)',(itype(i),i=1,nres)
294
295       do i=1,nres
296 #ifdef PROCOR
297         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
298 #else
299         if (itype(i).eq.21) then
300 #endif
301           itel(i)=0
302 #ifdef PROCOR
303         else if (itype(i+1).ne.20) then
304 #else
305         else if (itype(i).ne.20) then
306 #endif
307           itel(i)=1
308         else
309           itel(i)=2
310         endif
311       enddo
312       write (iout,*) "ITEL"
313       do i=1,nres-1
314         write (iout,*) i,itype(i),itel(i)
315       enddo
316
317       print *,'Call Read_Bridge.'
318       call read_bridge
319       if (with_dihed_constr) then
320
321       read (inp,*) ndih_constr
322       if (ndih_constr.gt.0) then
323         read (inp,*) ftors
324         write (iout,*) 'FTORS',ftors
325         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
326         write (iout,*)
327      &   'There are',ndih_constr,' constraints on phi angles.'
328         do i=1,ndih_constr
329           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
330         enddo
331         do i=1,ndih_constr
332           phi0(i)=deg2rad*phi0(i)
333           drange(i)=deg2rad*drange(i)
334         enddo
335       endif
336
337       endif
338
339       nnt=1
340       nct=nres
341       print *,'NNT=',NNT,' NCT=',NCT
342       if (itype(1).eq.21) nnt=2
343       if (itype(nres).eq.21) nct=nct-1
344       if (nstart.lt.nnt) nstart=nnt
345       if (nend.gt.nct .or. nend.eq.0) nend=nct
346       write (iout,*) "nstart",nstart," nend",nend
347       nres0=nres
348
349 C     Juyong:READ init_vars
350 C     Initialize variables!
351 C     Juyong:READ read_info
352 C     READ fragment information!!
353 C     both routines should be in dfa.F file!!
354
355       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
356      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
357 #ifdef DEBUG
358        write (iout,*) "Calling init_dfa_vars"
359 #ifdef AIX
360       call flush_(iout)
361 #else
362       call flush(iout)
363 #endif
364 #endif
365        call init_dfa_vars
366 #ifdef DEBUG
367        write (iout,*) 'init_dfa_vars finished!'
368 #ifdef AIX
369       call flush_(iout)
370 #else
371       call flush(iout)
372 #endif
373 #endif
374        call read_dfa_info
375 #ifdef DEBUG
376        write (iout,*) 'read_dfa_info finished!'
377 #ifdef AIX
378       call flush_(iout)
379 #else
380       call flush(iout)
381 #endif
382 #endif
383       endif
384
385       if (constr_homology.gt.0) then
386         call read_constr_homology(print_homology_restraints)
387       endif
388
389 c      if (pdbref) then
390 c        read(inp,'(a)') pdbfile
391 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
392 c        open(ipdbin,file=pdbfile,status='old',err=33)
393 c        goto 34 
394 c  33    write (iout,'(a)') 'Error opening PDB file.'
395 c        stop
396 c  34    continue
397 c        print *,'Begin reading pdb data'
398 c        call readpdb
399 c        print *,'Finished reading pdb data'
400 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
401 c        do i=1,nres
402 c          itype_pdb(i)=itype(i)
403 c        enddo
404 c        close (ipdbin)
405 c        write (iout,'(a,i3)') 'nsup=',nsup
406 c        nstart_seq=nnt
407 c        if (nsup.le.(nct-nnt+1)) then
408 c          do i=0,nct-nnt+1-nsup
409 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
410 c              nstart_seq=nnt+i
411 c              goto 111
412 c            endif
413 c          enddo
414 c          write (iout,'(a)') 
415 c     &            'Error - sequences to be superposed do not match.'
416 c          stop
417 c        else
418 c          do i=0,nsup-(nct-nnt+1)
419 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
420 c     &      then
421 c              nstart_sup=nstart_sup+i
422 c              nsup=nct-nnt+1
423 c              goto 111
424 c            endif
425 c          enddo 
426 c          write (iout,'(a)') 
427 c     &            'Error - sequences to be superposed do not match.'
428 c        endif
429 c  111   continue
430 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
431 c     &                 ' nstart_seq=',nstart_seq
432 c      endif
433       call init_int_table
434       call setup_var
435       write (iout,*) "molread: REFSTR",refstr
436       if (refstr) then
437         if (.not.pdbref) then
438           call read_angles(inp,*38)
439           goto 39
440    38     write (iout,'(a)') 'Error reading reference structure.'
441 #ifdef MPL
442           call mp_stopall(Error_Msg)
443 #else
444           stop 'Error reading reference structure'
445 #endif
446    39     call chainbuild     
447           nstart_sup=nnt
448           nstart_seq=nnt
449           nsup=nct-nnt+1
450           do i=1,2*nres
451             do j=1,3
452               cref(j,i)=c(j,i)
453             enddo
454           enddo
455         endif
456         call contact(print_contact_map,ncont_ref,icont_ref)
457       endif
458 c Read distance restraints
459       if (constr_dist.gt.0) then
460         call read_dist_constr
461         call hpb_partition
462       endif
463       return
464       end
465 c-----------------------------------------------------------------------------
466       logical function seq_comp(itypea,itypeb,length)
467       implicit none
468       integer length,itypea(length),itypeb(length)
469       integer i
470       do i=1,length
471         if (itypea(i).ne.itypeb(i)) then
472           seq_comp=.false.
473           return
474         endif
475       enddo
476       seq_comp=.true.
477       return
478       end
479 c-----------------------------------------------------------------------------
480       subroutine read_bridge
481 C Read information about disulfide bridges.
482       implicit none
483       include 'DIMENSIONS'
484       include 'COMMON.IOUNITS'
485       include 'COMMON.GEO'
486       include 'COMMON.VAR'
487       include 'COMMON.INTERACT'
488       include 'COMMON.LOCAL'
489       include 'COMMON.NAMES'
490       include 'COMMON.CHAIN'
491       include 'COMMON.FFIELD'
492       include 'COMMON.SBRIDGE'
493       include 'COMMON.HEADER'
494       include 'COMMON.CONTROL'
495       include 'COMMON.TIME1'
496 #ifdef MPL
497       include 'COMMON.INFO'
498 #endif
499       integer i,j
500 C Read bridging residues.
501       read (inp,*) ns,(iss(i),i=1,ns)
502       write(iout,*)'ns=',ns
503 C Check whether the specified bridging residues are cystines.
504       do i=1,ns
505         if (itype(iss(i)).ne.1) then
506           write (iout,'(2a,i3,a)') 
507      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
508      &   ' can form a disulfide bridge?!!!'
509           write (*,'(2a,i3,a)') 
510      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
511      &   ' can form a disulfide bridge?!!!'
512 #ifdef MPL
513          call mp_stopall(error_msg)
514 #else
515          stop
516 #endif
517         endif
518       enddo
519 C Read preformed bridges.
520       if (ns.gt.0) then
521       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
522       if (nss.gt.0) then
523         nhpb=nss
524 C Check if the residues involved in bridges are in the specified list of
525 C bridging residues.
526         do i=1,nss
527           do j=1,i-1
528             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
529      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
530               write (iout,'(a,i3,a)') 'Disulfide pair',i,
531      &      ' contains residues present in other pairs.'
532               write (*,'(a,i3,a)') 'Disulfide pair',i,
533      &      ' contains residues present in other pairs.'
534 #ifdef MPL
535               call mp_stopall(error_msg)
536 #else
537               stop 
538 #endif
539             endif
540           enddo
541           do j=1,ns
542             if (ihpb(i).eq.iss(j)) goto 10
543           enddo
544           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
545    10     continue
546           do j=1,ns
547             if (jhpb(i).eq.iss(j)) goto 20
548           enddo
549           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
550    20     continue
551 c          dhpb(i)=dbr
552 c          forcon(i)=fbr
553         enddo
554         do i=1,nss
555           ihpb(i)=ihpb(i)+nres
556           jhpb(i)=jhpb(i)+nres
557         enddo
558       endif
559       endif
560       if (ns.gt.0.and.dyn_ss) then
561           do i=nss+1,nhpb
562             ihpb(i-nss)=ihpb(i)
563             jhpb(i-nss)=jhpb(i)
564             forcon(i-nss)=forcon(i)
565             dhpb(i-nss)=dhpb(i)
566           enddo
567           nhpb=nhpb-nss
568           nss=0
569           call hpb_partition
570           do i=1,ns
571             dyn_ss_mask(iss(i))=.true.
572 c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
573           enddo
574       endif
575       print *, "Leaving brigde read"
576       return
577       end
578 c----------------------------------------------------------------------------
579       subroutine read_angles(kanal,*)
580       implicit none
581       include 'DIMENSIONS'
582       include 'COMMON.GEO'
583       include 'COMMON.VAR'
584       include 'COMMON.CHAIN'
585       include 'COMMON.IOUNITS'
586       integer i,kanal
587       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
588       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
589       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
590       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
591       do i=1,nres
592         theta(i)=deg2rad*theta(i)
593         phi(i)=deg2rad*phi(i)
594         alph(i)=deg2rad*alph(i)
595         omeg(i)=deg2rad*omeg(i)
596       enddo
597       return
598    10 return1
599       end
600 c----------------------------------------------------------------------------
601       subroutine reada(rekord,lancuch,wartosc,default)
602       implicit none
603       character*(*) rekord,lancuch
604       double precision wartosc,default
605       integer ilen,iread
606       external ilen
607       iread=index(rekord,lancuch)
608       if (iread.eq.0) then
609         wartosc=default 
610         return
611       endif   
612       iread=iread+ilen(lancuch)+1
613       read (rekord(iread:),*) wartosc
614       return
615       end
616 c----------------------------------------------------------------------------
617       subroutine multreada(rekord,lancuch,tablica,dim,default)
618       implicit none
619       integer dim,i
620       double precision tablica(dim),default
621       character*(*) rekord,lancuch
622       integer ilen,iread
623       external ilen
624       do i=1,dim
625         tablica(i)=default 
626       enddo
627       iread=index(rekord,lancuch)
628       if (iread.eq.0) return
629       iread=iread+ilen(lancuch)+1
630       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
631    10 return
632       end
633 c----------------------------------------------------------------------------
634       subroutine readi(rekord,lancuch,wartosc,default)
635       implicit none
636       character*(*) rekord,lancuch
637       integer wartosc,default
638       integer ilen,iread
639       external ilen
640       iread=index(rekord,lancuch)
641       if (iread.eq.0) then
642         wartosc=default 
643         return
644       endif   
645       iread=iread+ilen(lancuch)+1
646       read (rekord(iread:),*) wartosc
647       return
648       end
649 c----------------------------------------------------------------------------
650       subroutine multreadi(rekord,lancuch,tablica,dim,default)
651       implicit none
652       integer dim,i
653       integer tablica(dim),default
654       character*(*) rekord,lancuch
655       character*80 aux
656       integer ilen,iread
657       external ilen
658       do i=1,dim
659         tablica(i)=default
660       enddo
661       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
662       if (iread.eq.0) return
663       iread=iread+ilen(lancuch)+1
664       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
665    10 return
666       end
667 c----------------------------------------------------------------------------
668       subroutine card_concat(card)
669       include 'DIMENSIONS'
670       include 'COMMON.IOUNITS'
671       character*(*) card
672       character*80 karta,ucase
673       external ilen
674       read (inp,'(a)') karta
675       karta=ucase(karta)
676       card=' '
677       do while (karta(80:80).eq.'&')
678         card=card(:ilen(card)+1)//karta(:79)
679         read (inp,'(a)') karta
680         karta=ucase(karta)
681       enddo
682       card=card(:ilen(card)+1)//karta
683       return
684       end
685 c----------------------------------------------------------------------------
686       subroutine openunits
687       implicit none
688       include 'DIMENSIONS'    
689 #ifdef MPI
690       include "mpif.h"
691       character*3 liczba
692       include "COMMON.MPI"
693 #endif
694       include 'COMMON.IOUNITS'
695       include 'COMMON.CONTROL'
696       integer lenpre,lenpot,ilen
697       external ilen
698       character*16 cformat,cprint
699       character*16 ucase
700       integer lenint,lenout
701       call getenv('INPUT',prefix)
702       call getenv('OUTPUT',prefout)
703       call getenv('INTIN',prefintin)
704       call getenv('COORD',cformat)
705       call getenv('PRINTCOOR',cprint)
706       call getenv('SCRATCHDIR',scratchdir)
707       from_bx=.true.
708       from_cx=.false.
709       if (index(ucase(cformat),'CX').gt.0) then
710         from_cx=.true.
711         from_bx=.false.
712       endif
713       from_cart=.true.
714       lenpre=ilen(prefix)
715       lenout=ilen(prefout)
716       lenint=ilen(prefintin)
717 C Get the names and open the input files
718       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
719 #ifdef MPI
720       write (liczba,'(bz,i3.3)') me
721       outname=prefout(:lenout)//'_clust.out_'//liczba
722 #else
723       outname=prefout(:lenout)//'_clust.out'
724 #endif
725       if (from_bx) then
726         intinname=prefintin(:lenint)//'.bx'
727       else if (from_cx) then
728         intinname=prefintin(:lenint)//'.cx'
729       else
730         intinname=prefintin(:lenint)//'.int'
731       endif
732       rmsname=prefintin(:lenint)//'.rms'
733       open (jplot,file=prefout(:ilen(prefout))//'.tex',
734      &       status='unknown')
735       open (jrms,file=rmsname,status='unknown')
736       open(iout,file=outname,status='unknown')
737 C Get parameter filenames and open the parameter files.
738       call getenv('BONDPAR',bondname)
739       open (ibond,file=bondname,status='old')
740       call getenv('THETPAR',thetname)
741       open (ithep,file=thetname,status='old')
742       call getenv('ROTPAR',rotname)
743       open (irotam,file=rotname,status='old')
744       call getenv('TORPAR',torname)
745       open (itorp,file=torname,status='old')
746       call getenv('TORDPAR',tordname)
747       open (itordp,file=tordname,status='old')
748       call getenv('FOURIER',fouriername)
749       open (ifourier,file=fouriername,status='old')
750       call getenv('ELEPAR',elename)
751       open (ielep,file=elename,status='old')
752       call getenv('SIDEPAR',sidename)
753       open (isidep,file=sidename,status='old')
754       call getenv('SIDEP',sidepname)
755       open (isidep1,file=sidepname,status="old")
756       call getenv('SCCORPAR',sccorname)
757       open (isccor,file=sccorname,status="old")
758 #ifndef OLDSCP
759 C
760 C 8/9/01 In the newest version SCp interaction constants are read from a file
761 C Use -DOLDSCP to use hard-coded constants instead.
762 C
763       call getenv('SCPPAR',scpname)
764       open (iscpp,file=scpname,status='old')
765 #endif
766       return
767       end
768 c-------------------------------------------------------------------------------
769       subroutine read_dist_constr
770       implicit real*8 (a-h,o-z)
771       include 'DIMENSIONS'
772 #ifdef MPI
773       include 'mpif.h'
774 #endif
775       include 'COMMON.CONTROL'
776       include 'COMMON.CHAIN'
777       include 'COMMON.IOUNITS'
778       include 'COMMON.SBRIDGE'
779       integer ifrag_(2,100),ipair_(2,100)
780       double precision wfrag_(100),wpair_(100)
781       character*500 controlcard
782       logical lprn /.true./
783       logical normalize,next
784       integer restr_type
785       double precision xlink(4,0:4) /
786 c           a          b       c     sigma
787      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
788      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
789      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
790      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
791      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
792       write (iout,*) "Calling read_dist_constr"
793 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
794 c      call flush(iout)
795       next=.true.
796
797       DO WHILE (next)
798
799       call card_concat(controlcard)
800       next = index(controlcard,"NEXT").gt.0
801       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
802       write (iout,*) "restr_type",restr_type
803       call readi(controlcard,"NFRAG",nfrag_,0)
804       call readi(controlcard,"NFRAG",nfrag_,0)
805       call readi(controlcard,"NPAIR",npair_,0)
806       call readi(controlcard,"NDIST",ndist_,0)
807       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
808       if (restr_type.eq.10) 
809      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
810       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
811       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
812       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
813       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
814       normalize = index(controlcard,"NORMALIZE").gt.0
815       write (iout,*) "WBOLTZD",wboltzd
816 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
817 c      write (iout,*) "IFRAG"
818 c      do i=1,nfrag_
819 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
820 c      enddo
821 c      write (iout,*) "IPAIR"
822 c      do i=1,npair_
823 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
824 c      enddo
825       if (nfrag_.gt.0) then
826         nres0=nres
827         read(inp,'(a)') pdbfile
828         write (iout,*) 
829      & "Distance restraints will be constructed from structure ",pdbfile
830         open(ipdbin,file=pdbfile,status='old',err=11)
831         call readpdb(.true.)
832         nres=nres0
833         close(ipdbin)
834       endif
835       do i=1,nfrag_
836         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
837         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
838      &    ifrag_(2,i)=nstart_sup+nsup-1
839 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
840 c        call flush(iout)
841         if (wfrag_(i).eq.0.0d0) cycle
842         do j=ifrag_(1,i),ifrag_(2,i)-1
843           do k=j+1,ifrag_(2,i)
844 c            write (iout,*) "j",j," k",k
845             ddjk=dist(j,k)
846             if (restr_type.eq.1) then
847               nhpb=nhpb+1
848               irestr_type(nhpb)=1
849               ihpb(nhpb)=j
850               jhpb(nhpb)=k
851               dhpb(nhpb)=ddjk
852               forcon(nhpb)=wfrag_(i) 
853             else if (constr_dist.eq.2) then
854               if (ddjk.le.dist_cut) then
855                 nhpb=nhpb+1
856                 irestr_type(nhpb)=1
857                 ihpb(nhpb)=j
858                 jhpb(nhpb)=k
859                 dhpb(nhpb)=ddjk
860                 forcon(nhpb)=wfrag_(i) 
861               endif
862             else if (restr_type.eq.3) then
863               nhpb=nhpb+1
864               irestr_type(nhpb)=1
865               ihpb(nhpb)=j
866               jhpb(nhpb)=k
867               dhpb(nhpb)=ddjk
868               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
869             endif
870             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
871      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
872           enddo
873         enddo
874       enddo
875       do i=1,npair_
876         if (wpair_(i).eq.0.0d0) cycle
877         ii = ipair_(1,i)
878         jj = ipair_(2,i)
879         if (ii.gt.jj) then
880           itemp=ii
881           ii=jj
882           jj=itemp
883         endif
884         do j=ifrag_(1,ii),ifrag_(2,ii)
885           do k=ifrag_(1,jj),ifrag_(2,jj)
886             if (restr_type.eq.1) then
887               nhpb=nhpb+1
888               irestr_type(nhpb)=1
889               ihpb(nhpb)=j
890               jhpb(nhpb)=k
891               dhpb(nhpb)=ddjk
892               forcon(nhpb)=wfrag_(i) 
893             else if (constr_dist.eq.2) then
894               if (ddjk.le.dist_cut) then
895                 nhpb=nhpb+1
896                 irestr_type(nhpb)=1
897                 ihpb(nhpb)=j
898                 jhpb(nhpb)=k
899                 dhpb(nhpb)=ddjk
900                 forcon(nhpb)=wfrag_(i) 
901               endif
902             else if (restr_type.eq.3) then
903               nhpb=nhpb+1
904               irestr_type(nhpb)=1
905               ihpb(nhpb)=j
906               jhpb(nhpb)=k
907               dhpb(nhpb)=ddjk
908               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
909             endif
910             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
911      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
912           enddo
913         enddo
914       enddo 
915
916 c      print *,ndist_
917       write (iout,*) "Distance restraints as read from input"
918       do i=1,ndist_
919         if (restr_type.eq.11) then
920           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
921      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
922 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
923           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
924           nhpb=nhpb+1
925           irestr_type(nhpb)=11
926           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
927      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
928      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
929           if (ibecarb(nhpb).gt.0) then
930             ihpb(nhpb)=ihpb(nhpb)+nres
931             jhpb(nhpb)=jhpb(nhpb)+nres
932           endif
933         else if (constr_dist.eq.10) then
934 c Cross-lonk Markov-like potential
935           call card_concat(controlcard)
936           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
937           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
938           ibecarb(nhpb+1)=0
939           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
940           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
941           if (index(controlcard,"ZL").gt.0) then
942             link_type=1
943           else if (index(controlcard,"ADH").gt.0) then
944             link_type=2
945           else if (index(controlcard,"PDH").gt.0) then
946             link_type=3
947           else if (index(controlcard,"DSS").gt.0) then
948             link_type=4
949           else
950             link_type=0
951           endif
952           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
953      &       xlink(1,link_type))
954           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
955      &       xlink(2,link_type))
956           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
957      &       xlink(3,link_type))
958           call reada(controlcard,"SIGMA",forcon(nhpb+1),
959      &       xlink(4,link_type))
960           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
961 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
962 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
963           if (forcon(nhpb+1).le.0.0d0 .or. 
964      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
965           nhpb=nhpb+1
966           irestr_type(nhpb)=10
967           if (ibecarb(nhpb).gt.0) then
968             ihpb(nhpb)=ihpb(nhpb)+nres
969             jhpb(nhpb)=jhpb(nhpb)+nres
970           endif
971           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
972      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
973      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
974      &     irestr_type(nhpb)
975         else
976 C        print *,"in else"
977           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
978      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
979           if (forcon(nhpb+1).gt.0.0d0) then
980           nhpb=nhpb+1
981           if (dhpb1(nhpb).eq.0.0d0) then
982             irestr_type(nhpb)=1
983           else
984             irestr_type(nhpb)=2
985           endif
986           if (ibecarb(nhpb).gt.0) then
987             ihpb(nhpb)=ihpb(nhpb)+nres
988             jhpb(nhpb)=jhpb(nhpb)+nres
989           endif
990           if (dhpb(nhpb).eq.0.0d0)
991      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
992         endif
993         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
994      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
995         endif
996 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
997 C        if (forcon(nhpb+1).gt.0.0d0) then
998 C          nhpb=nhpb+1
999 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1000       enddo
1001
1002       ENDDO ! next
1003
1004       fordepthmax=0.0d0
1005       if (normalize) then
1006         do i=nss+1,nhpb
1007           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax) 
1008      &      fordepthmax=fordepth(i)
1009         enddo
1010         do i=nss+1,nhpb
1011           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1012         enddo
1013       endif
1014       if (nhpb.gt.nss)  then
1015         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1016      &  "The following",nhpb-nss,
1017      &  " distance restraints have been imposed:",
1018      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
1019      &  "  score"," type"
1020         do i=nss+1,nhpb
1021           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1022      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1023      &  irestr_type(i)
1024         enddo
1025       endif
1026       call hpb_partition
1027       call flush(iout)
1028       return
1029    11 write (iout,*)"read_dist_restr: error reading reference structure"
1030       stop
1031       end
1032 c====-------------------------------------------------------------------
1033       subroutine read_constr_homology(lprn)
1034
1035       include 'DIMENSIONS'
1036 #ifdef MPI
1037       include 'mpif.h'
1038 #endif
1039       include 'COMMON.SETUP'
1040       include 'COMMON.CONTROL'
1041       include 'COMMON.CHAIN'
1042       include 'COMMON.IOUNITS'
1043       include 'COMMON.GEO'
1044       include 'COMMON.INTERACT'
1045       include 'COMMON.NAMES'
1046       include 'COMMON.HOMRESTR'
1047 c
1048 c For new homol impl
1049 c
1050       include 'COMMON.VAR'
1051 c     include 'include_unres/COMMON.VAR'
1052 c
1053
1054 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1055 c    &                 dist_cut
1056 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1057 c    &    sigma_odl_temp(maxres,maxres,max_template)
1058       character*2 kic2
1059       character*24 model_ki_dist, model_ki_angle
1060       character*500 controlcard
1061       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1062       integer idomain(max_template,maxres)
1063       integer ilen
1064       external ilen
1065       logical lprn
1066       logical unres_pdb,liiflag
1067 c
1068 c     FP - Nov. 2014 Temporary specifications for new vars
1069 c
1070       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1071       double precision, dimension (max_template,maxres) :: rescore
1072       double precision, dimension (max_template,maxres) :: rescore2
1073       character*24 tpl_k_rescore
1074 c -----------------------------------------------------------------
1075 c Reading multiple PDB ref structures and calculation of retraints
1076 c not using pre-computed ones stored in files model_ki_{dist,angle}
1077 c FP (Nov., 2014)
1078 c -----------------------------------------------------------------
1079 c
1080 c
1081 c Alternative: reading from input
1082 #ifdef DEBUG
1083       write (iout,*) "BEGIN READ HOMOLOGY INFO"
1084 #ifdef AIX
1085       call flush_(iout)
1086 #else
1087       call flush(iout)
1088 #endif
1089 #endif
1090       call card_concat(controlcard)
1091       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1092       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1093       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1094       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1095       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1096       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1097       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1098       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1099       if (homol_nset.gt.1)then
1100          call readi(controlcard,"ISET",iset,1)
1101          call card_concat(controlcard)
1102          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1103       else
1104         iset=1
1105         waga_homology(1)=1.0
1106       endif
1107 c
1108 #ifdef DEBUG
1109       write(iout,*) "read_constr_homology iset",iset
1110       write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
1111 #ifdef AIX
1112       call flush_(iout)
1113 #else
1114       call flush(iout)
1115 #endif
1116 #endif
1117 cd      write (iout,*) "nnt",nnt," nct",nct
1118 cd      call flush(iout)
1119
1120
1121       lim_odl=0
1122       lim_dih=0
1123 c
1124 c  Reading HM global scores (prob not required)
1125 c
1126       do i = nnt,nct
1127         do k=1,constr_homology
1128          idomain(k,i)=0
1129         enddo
1130       enddo
1131 c     open (4,file="HMscore")
1132 c     do k=1,constr_homology
1133 c       read (4,*,end=521) hmscore_tmp
1134 c       hmscore(k)=hmscore_tmp ! Another transformation can be used 
1135 c       write(*,*) "Model", k, ":", hmscore(k)
1136 c     enddo
1137 c521  continue
1138
1139       ii=0
1140       do i = nnt,nct-2 
1141         do j=i+2,nct 
1142         ii=ii+1
1143         ii_in_use(ii)=0
1144         enddo
1145       enddo
1146 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1147
1148       if (read_homol_frag) then   
1149         call read_klapaucjusz   
1150       else 
1151       write (iout,*) "CONSTR_HOMOLOGY",constr_homology
1152       do k=1,constr_homology
1153
1154         read(inp,'(a)') pdbfile
1155 c        write (iout,*) "k ",k," pdbfile ",pdbfile
1156 c  Next stament causes error upon compilation (?)
1157 c       if(me.eq.king.or. .not. out1file)
1158 c         write (iout,'(2a)') 'PDB data will be read from file ',
1159 c    &   pdbfile(:ilen(pdbfile))
1160          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
1161      &  pdbfile(:ilen(pdbfile))
1162         open(ipdbin,file=pdbfile,status='old',err=33)
1163         goto 34
1164   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1165      &  pdbfile(:ilen(pdbfile))
1166         stop
1167   34    continue
1168 c        print *,'Begin reading pdb data'
1169 c
1170 c Files containing res sim or local scores (former containing sigmas)
1171 c
1172
1173         write(kic2,'(bz,i2.2)') k
1174
1175         tpl_k_rescore="template"//kic2//".sco"
1176
1177         unres_pdb=.false.
1178         call readpdb
1179         do i=1,2*nres
1180           do j=1,3
1181             crefjlee(j,i)=c(j,i)
1182           enddo
1183         enddo
1184 #ifdef DEBUG
1185         do i=1,nres
1186           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1187      &      (crefjlee(j,i+nres),j=1,3)
1188         enddo
1189         write (iout,*) "READ HOMOLOGY INFO"
1190         write (iout,*) "read_constr_homology x: after reading pdb file"
1191         write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
1192         write (iout,*) "waga_dist",waga_dist
1193         write (iout,*) "waga_angle",waga_angle
1194         write (iout,*) "waga_theta",waga_theta
1195         write (iout,*) "waga_d",waga_d
1196         write (iout,*) "dist_cut",dist_cut
1197 #endif
1198 #ifdef AIX
1199       call flush_(iout)
1200 #else
1201       call flush(iout)
1202 #endif
1203
1204 c
1205 c     Distance restraints
1206 c
1207 c          ... --> odl(k,ii)
1208 C Copy the coordinates from reference coordinates (?)
1209         do i=1,2*nres
1210           do j=1,3
1211             c(j,i)=cref(j,i)
1212 c           write (iout,*) "c(",j,i,") =",c(j,i)
1213           enddo
1214         enddo
1215 c
1216 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
1217 c
1218 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1219           open (ientin,file=tpl_k_rescore,status='old')
1220           if (nnt.gt.1) rescore(k,1)=0.0d0
1221           do irec=nnt,nct ! loop for reading res sim 
1222             if (read2sigma) then
1223              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
1224      &                                idomain_tmp
1225              i_tmp=i_tmp+nnt-1
1226              idomain(k,i_tmp)=idomain_tmp
1227              rescore(k,i_tmp)=rescore_tmp
1228              rescore2(k,i_tmp)=rescore2_tmp
1229              write(iout,'(a7,i5,2f10.5,i5)') "rescore",
1230      &                      i_tmp,rescore2_tmp,rescore_tmp,
1231      &                                idomain_tmp
1232             else
1233              idomain(k,irec)=1
1234              read (ientin,*,end=1401) rescore_tmp
1235
1236 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1237              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1238 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1239             endif
1240           enddo  
1241  1401   continue
1242         close (ientin)        
1243         if (waga_dist.ne.0.0d0) then
1244           ii=0
1245           do i = nnt,nct-2 
1246             do j=i+2,nct 
1247
1248               x12=c(1,i)-c(1,j)
1249               y12=c(2,i)-c(2,j)
1250               z12=c(3,i)-c(3,j)
1251               distal=dsqrt(x12*x12+y12*y12+z12*z12) 
1252 c              write (iout,*) k,i,j,distal,dist2_cut
1253
1254             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1255      &            .and. distal.le.dist2_cut ) then
1256
1257               ii=ii+1
1258               ii_in_use(ii)=1
1259               l_homo(k,ii)=.true.
1260
1261 c             write (iout,*) "k",k
1262 c             write (iout,*) "i",i," j",j," constr_homology",
1263 c    &                       constr_homology
1264               ires_homo(ii)=i
1265               jres_homo(ii)=j
1266               odl(k,ii)=distal
1267               if (read2sigma) then
1268                 sigma_odl(k,ii)=0
1269                 do ik=i,j
1270                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1271                 enddo
1272                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1273                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = 
1274      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1275               else
1276                 if (odl(k,ii).le.dist_cut) then
1277                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) 
1278                 else
1279 #ifdef OLDSIGMA
1280                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
1281      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1282 #else
1283                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
1284      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1285 #endif
1286                 endif
1287               endif
1288               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) 
1289             else
1290               ii=ii+1
1291               l_homo(k,ii)=.false.
1292             endif
1293             enddo
1294           enddo
1295         lim_odl=ii
1296         endif
1297 c
1298 c     Theta, dihedral and SC retraints
1299 c
1300         if (waga_angle.gt.0.0d0) then
1301 c         open (ientin,file=tpl_k_sigma_dih,status='old')
1302 c         do irec=1,maxres-3 ! loop for reading sigma_dih
1303 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1304 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1305 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1306 c    &                            sigma_dih(k,i+nnt-1)
1307 c         enddo
1308 c1402   continue
1309 c         close (ientin)
1310           do i = nnt+3,nct 
1311             if (idomain(k,i).eq.0) then 
1312                sigma_dih(k,i)=0.0
1313                cycle
1314             endif
1315             dih(k,i)=phiref(i) ! right?
1316 c           read (ientin,*) sigma_dih(k,i) ! original variant
1317 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
1318 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1319 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1320 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
1321 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
1322
1323             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1324      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1325 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1326 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
1327 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1328 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1329 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
1330 c   Instead of res sim other local measure of b/b str reliability possible
1331             if (sigma_dih(k,i).ne.0)
1332      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1333 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1334           enddo
1335           lim_dih=nct-nnt-2 
1336         endif
1337
1338         if (waga_theta.gt.0.0d0) then
1339 c         open (ientin,file=tpl_k_sigma_theta,status='old')
1340 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1341 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1342 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1343 c    &                              sigma_theta(k,i+nnt-1)
1344 c         enddo
1345 c1403   continue
1346 c         close (ientin)
1347
1348           do i = nnt+2,nct ! right? without parallel.
1349 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
1350 c         do i=ithet_start,ithet_end ! with FG parallel.
1351              if (idomain(k,i).eq.0) then  
1352               sigma_theta(k,i)=0.0
1353               cycle
1354              endif
1355              thetatpl(k,i)=thetaref(i)
1356 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1357 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
1358 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
1359 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
1360 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
1361              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1362      &                        rescore(k,i-2))/3.0
1363 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1364              if (sigma_theta(k,i).ne.0)
1365      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1366
1367 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1368 c                             rescore(k,i-2) !  right expression ?
1369 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1370           enddo
1371         endif
1372
1373         if (waga_d.gt.0.0d0) then
1374 c       open (ientin,file=tpl_k_sigma_d,status='old')
1375 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1376 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1377 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1378 c    &                          sigma_d(k,i+nnt-1)
1379 c         enddo
1380 c1404   continue
1381
1382           do i = nnt,nct ! right? without parallel.
1383 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
1384 c         do i=loc_start,loc_end ! with FG parallel.
1385                if (itype(i).eq.10) cycle 
1386                if (idomain(k,i).eq.0 ) then 
1387                   sigma_d(k,i)=0.0
1388                   cycle
1389                endif
1390                xxtpl(k,i)=xxref(i)
1391                yytpl(k,i)=yyref(i)
1392                zztpl(k,i)=zzref(i)
1393 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1394 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1395 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1396 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
1397                sigma_d(k,i)=rescore(k,i) !  right expression ?
1398                if (sigma_d(k,i).ne.0)
1399      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1400
1401 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
1402 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1403 c              read (ientin,*) sigma_d(k,i) ! 1st variant
1404                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
1405           enddo
1406         endif
1407       enddo
1408 c
1409 c remove distance restraints not used in any model from the list
1410 c shift data in all arrays
1411 c
1412       if (waga_dist.ne.0.0d0) then
1413         ii=0
1414         liiflag=.true.
1415         do i=nnt,nct-2 
1416          do j=i+2,nct 
1417           ii=ii+1
1418           if (ii_in_use(ii).eq.0.and.liiflag) then
1419             liiflag=.false.
1420             iistart=ii
1421           endif
1422           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1423      &                   .not.liiflag.and.ii.eq.lim_odl) then
1424              if (ii.eq.lim_odl) then
1425               iishift=ii-iistart+1
1426              else
1427               iishift=ii-iistart
1428              endif
1429              liiflag=.true.
1430              do ki=iistart,lim_odl-iishift
1431               ires_homo(ki)=ires_homo(ki+iishift)
1432               jres_homo(ki)=jres_homo(ki+iishift)
1433               ii_in_use(ki)=ii_in_use(ki+iishift)
1434               do k=1,constr_homology
1435                odl(k,ki)=odl(k,ki+iishift)
1436                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1437                l_homo(k,ki)=l_homo(k,ki+iishift)
1438               enddo
1439              enddo
1440              ii=ii-iishift
1441              lim_odl=lim_odl-iishift
1442           endif
1443          enddo
1444         enddo
1445       endif
1446
1447       endif ! .not. klapaucjusz
1448
1449       if (constr_homology.gt.0) call homology_partition
1450       if (constr_homology.gt.0) call init_int_table
1451 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
1452 cd     & "lim_xx=",lim_xx
1453 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1454 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1455 c
1456 c Print restraints
1457 c
1458       if (.not.lprn) return
1459 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1460       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1461        write (iout,*) "Distance restraints from templates"
1462        do ii=1,lim_odl
1463        write(iout,'(3i5,100(2f8.2,1x,l1,4x))') 
1464      &  ii,ires_homo(ii),jres_homo(ii),
1465      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
1466      &  ki=1,constr_homology)
1467        enddo
1468        write (iout,*) "Dihedral angle restraints from templates"
1469        do i=nnt+3,nct
1470         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1471      &      (rad2deg*dih(ki,i),
1472      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1473        enddo
1474        write (iout,*) "Virtual-bond angle restraints from templates"
1475        do i=nnt+2,nct
1476         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
1477      &      (rad2deg*thetatpl(ki,i),
1478      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1479        enddo
1480        write (iout,*) "SC restraints from templates"
1481        do i=nnt,nct
1482         write(iout,'(i5,100(4f8.2,4x))') i,
1483      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
1484      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1485        enddo
1486       endif
1487 c -----------------------------------------------------------------
1488       return
1489       end
1490 c----------------------------------------------------------------------
1491       subroutine read_klapaucjusz
1492
1493       include 'DIMENSIONS'
1494 #ifdef MPI
1495       include 'mpif.h'
1496 #endif
1497       include 'COMMON.SETUP'
1498       include 'COMMON.CONTROL'
1499       include 'COMMON.CHAIN'
1500       include 'COMMON.IOUNITS'
1501       include 'COMMON.GEO'
1502       include 'COMMON.INTERACT'
1503       include 'COMMON.NAMES'
1504       include 'COMMON.HOMRESTR'
1505       character*256 fragfile
1506       integer ninclust(maxclust),inclust(max_template,maxclust),
1507      &  nresclust(maxclust),iresclust(maxres,maxclust)
1508
1509       character*2 kic2
1510       character*24 model_ki_dist, model_ki_angle
1511       character*500 controlcard
1512       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
1513       integer idomain(max_template,maxres)
1514       logical lprn /.true./
1515       integer ilen
1516       external ilen
1517       logical unres_pdb,liiflag
1518 c
1519 c
1520       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
1521       double precision, dimension (max_template,maxres) :: rescore
1522       double precision, dimension (max_template,maxres) :: rescore2
1523       character*24 tpl_k_rescore
1524
1525 c
1526 c For new homol impl
1527 c
1528       include 'COMMON.VAR'
1529 c
1530       double precision chomo(3,maxres2+2,max_template)
1531       call getenv("FRAGFILE",fragfile) 
1532       open(ientin,file=fragfile,status="old",err=10)
1533       read(ientin,*) constr_homology,nclust
1534       l_homo = .false.
1535       sigma_theta=0.0
1536       sigma_d=0.0
1537       sigma_dih=0.0
1538 c Read pdb files 
1539       do k=1,constr_homology 
1540         read(ientin,'(a)') pdbfile
1541         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
1542      &  pdbfile(:ilen(pdbfile))
1543         open(ipdbin,file=pdbfile,status='old',err=33)
1544         goto 34
1545   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
1546      &  pdbfile(:ilen(pdbfile))
1547         stop
1548   34    continue
1549         unres_pdb=.false.
1550         call readpdb
1551         do i=1,2*nres
1552           do j=1,3
1553             chomo(j,i,k)=c(j,i)
1554           enddo
1555         enddo
1556         do i=1,nres
1557           rescore(k,i)=0.2d0
1558           rescore2(k,i)=1.0d0
1559         enddo
1560       enddo
1561 c Read clusters
1562       do i=1,nclust
1563         read(ientin,*) ninclust(i),nresclust(i)
1564         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
1565         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
1566       enddo
1567 c
1568 c Loop over clusters
1569 c
1570       do l=1,nclust
1571         do ll = 1,ninclust(l)
1572         
1573         k = inclust(ll,l)
1574         do i=1,nres
1575           idomain(k,i)=0
1576         enddo
1577         do i=1,nresclust(l)
1578           if (nnt.gt.1)  then
1579             idomain(k,iresclust(i,l)+1) = 1
1580           else
1581             idomain(k,iresclust(i,l)) = 1
1582           endif
1583         enddo
1584 c
1585 c     Distance restraints
1586 c
1587 c          ... --> odl(k,ii)
1588 C Copy the coordinates from reference coordinates (?)
1589         do i=1,2*nres
1590           do j=1,3
1591             c(j,i)=chomo(j,i,k)
1592 c           write (iout,*) "c(",j,i,") =",c(j,i)
1593           enddo
1594         enddo
1595         call int_from_cart(.true.,.false.)
1596         call sc_loc_geom(.false.)
1597         do i=1,nres
1598           thetaref(i)=theta(i)
1599           phiref(i)=phi(i)
1600         enddo
1601         if (waga_dist.ne.0.0d0) then
1602           ii=0
1603           do i = nnt,nct-2 
1604             do j=i+2,nct 
1605
1606               x12=c(1,i)-c(1,j)
1607               y12=c(2,i)-c(2,j)
1608               z12=c(3,i)-c(3,j)
1609               distal=dsqrt(x12*x12+y12*y12+z12*z12) 
1610 c              write (iout,*) k,i,j,distal,dist2_cut
1611
1612             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1613      &            .and. distal.le.dist2_cut ) then
1614
1615               ii=ii+1
1616               ii_in_use(ii)=1
1617               l_homo(k,ii)=.true.
1618
1619 c             write (iout,*) "k",k
1620 c             write (iout,*) "i",i," j",j," constr_homology",
1621 c    &                       constr_homology
1622               ires_homo(ii)=i
1623               jres_homo(ii)=j
1624               odl(k,ii)=distal
1625               if (read2sigma) then
1626                 sigma_odl(k,ii)=0
1627                 do ik=i,j
1628                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1629                 enddo
1630                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1631                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = 
1632      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1633               else
1634                 if (odl(k,ii).le.dist_cut) then
1635                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) 
1636                 else
1637 #ifdef OLDSIGMA
1638                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
1639      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1640 #else
1641                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* 
1642      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1643 #endif
1644                 endif
1645               endif
1646               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) 
1647             else
1648               ii=ii+1
1649 c              l_homo(k,ii)=.false.
1650             endif
1651             enddo
1652           enddo
1653         lim_odl=ii
1654         endif
1655 c
1656 c     Theta, dihedral and SC retraints
1657 c
1658         if (waga_angle.gt.0.0d0) then
1659           do i = nnt+3,nct 
1660             if (idomain(k,i).eq.0) then 
1661 c               sigma_dih(k,i)=0.0
1662                cycle
1663             endif
1664             dih(k,i)=phiref(i)
1665             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1666      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
1667 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1668 c     &       " sigma_dihed",sigma_dih(k,i)
1669             if (sigma_dih(k,i).ne.0)
1670      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1671           enddo
1672           lim_dih=nct-nnt-2 
1673         endif
1674
1675         if (waga_theta.gt.0.0d0) then
1676           do i = nnt+2,nct
1677              if (idomain(k,i).eq.0) then  
1678 c              sigma_theta(k,i)=0.0
1679               cycle
1680              endif
1681              thetatpl(k,i)=thetaref(i)
1682              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1683      &                        rescore(k,i-2))/3.0
1684              if (sigma_theta(k,i).ne.0)
1685      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1686           enddo
1687         endif
1688
1689         if (waga_d.gt.0.0d0) then
1690           do i = nnt,nct
1691                if (itype(i).eq.10) cycle 
1692                if (idomain(k,i).eq.0 ) then 
1693 c                  sigma_d(k,i)=0.0
1694                   cycle
1695                endif
1696                xxtpl(k,i)=xxref(i)
1697                yytpl(k,i)=yyref(i)
1698                zztpl(k,i)=zzref(i)
1699                sigma_d(k,i)=rescore(k,i)
1700                if (sigma_d(k,i).ne.0)
1701      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1702                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1703           enddo
1704         endif
1705       enddo ! l
1706       enddo ! ll
1707 c
1708 c remove distance restraints not used in any model from the list
1709 c shift data in all arrays
1710 c
1711       if (waga_dist.ne.0.0d0) then
1712         ii=0
1713         liiflag=.true.
1714         do i=nnt,nct-2 
1715          do j=i+2,nct 
1716           ii=ii+1
1717           if (ii_in_use(ii).eq.0.and.liiflag) then
1718             liiflag=.false.
1719             iistart=ii
1720           endif
1721           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1722      &                   .not.liiflag.and.ii.eq.lim_odl) then
1723              if (ii.eq.lim_odl) then
1724               iishift=ii-iistart+1
1725              else
1726               iishift=ii-iistart
1727              endif
1728              liiflag=.true.
1729              do ki=iistart,lim_odl-iishift
1730               ires_homo(ki)=ires_homo(ki+iishift)
1731               jres_homo(ki)=jres_homo(ki+iishift)
1732               ii_in_use(ki)=ii_in_use(ki+iishift)
1733               do k=1,constr_homology
1734                odl(k,ki)=odl(k,ki+iishift)
1735                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1736                l_homo(k,ki)=l_homo(k,ki+iishift)
1737               enddo
1738              enddo
1739              ii=ii-iishift
1740              lim_odl=lim_odl-iishift
1741           endif
1742          enddo
1743         enddo
1744       endif
1745
1746       return
1747    10 stop "Error infragment file"
1748       end
1749 c----------------------------------------------------------------------