Merge branch 'devel' of mmka:unres into devel
[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       call readi(controlcard,'RESCALE',rescale_mode,2)
28       call readi(controlcard,'PDBOUT',outpdb,0)
29       call readi(controlcard,'MOL2OUT',outmol2,0)
30       refstr=(index(controlcard,'REFSTR').gt.0)
31       write (iout,*) "REFSTR",refstr
32       pdbref=(index(controlcard,'PDBREF').gt.0)
33       iscode=index(controlcard,'ONE_LETTER')
34       tree=(index(controlcard,'MAKE_TREE').gt.0)
35       min_var=(index(controlcard,'MINVAR').gt.0)
36       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
37       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
38       call readi(controlcard,'NCUT',ncut,1)
39       call readi(controlcard,'NSTART',nstart,0)
40       call readi(controlcard,'NEND',nend,0)
41       call reada(controlcard,'ECUT',ecut,10.0d0)
42       call reada(controlcard,'PROB',prob_limit,0.99d0)
43       write (iout,*) "Probability limit",prob_limit
44       lgrp=(index(controlcard,'LGRP').gt.0)
45       caonly=(index(controlcard,'CA_ONLY').gt.0)
46       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
47       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
48       call readi(controlcard,'IOPT',iopt,2) 
49       lside = index(controlcard,"SIDE").gt.0
50       efree = index(controlcard,"EFREE").gt.0
51       call readi(controlcard,'NTEMP',nT,1)
52       write (iout,*) "nT",nT
53       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
54       write (iout,*) "nT",nT
55       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
56       do i=1,nT
57         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
58       enddo
59       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
60       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
61       lprint_int=index(controlcard,"PRINT_INT") .gt.0
62       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
63       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
64       write (iout,*) "with_dihed_constr ",with_dihed_constr,
65      & " CONSTR_DIST",constr_dist
66       call flush(iout)
67       if (min_var) iopt=1
68       return
69       end
70 c--------------------------------------------------------------------------
71       subroutine molread
72 C
73 C Read molecular data.
74 C
75       implicit none
76       include 'DIMENSIONS'
77       include 'COMMON.IOUNITS'
78       include 'COMMON.GEO'
79       include 'COMMON.VAR'
80       include 'COMMON.INTERACT'
81       include 'COMMON.LOCAL'
82       include 'COMMON.NAMES'
83       include 'COMMON.CHAIN'
84       include 'COMMON.FFIELD'
85       include 'COMMON.SBRIDGE'
86       include 'COMMON.HEADER'
87       include 'COMMON.CONTROL'
88       include 'COMMON.CONTACTS'
89       include 'COMMON.TIME1'
90       include 'COMMON.TORCNSTR'
91 #ifdef MPL
92       include 'COMMON.INFO'
93 #endif
94       character*4 sequence(maxres)
95       character*800 weightcard
96       integer rescode
97       double precision x(maxvar)
98       integer itype_pdb(maxres)
99       logical seq_comp
100       integer i,j
101 C
102 C Body
103 C
104 C Read weights of the subsequent energy terms.
105       call card_concat(weightcard)
106       call reada(weightcard,'WSC',wsc,1.0d0)
107       call reada(weightcard,'WLONG',wsc,wsc)
108       call reada(weightcard,'WSCP',wscp,1.0d0)
109       call reada(weightcard,'WELEC',welec,1.0D0)
110       call reada(weightcard,'WVDWPP',wvdwpp,welec)
111       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
112       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
113       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
114       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
115       call reada(weightcard,'WTURN3',wturn3,1.0D0)
116       call reada(weightcard,'WTURN4',wturn4,1.0D0)
117       call reada(weightcard,'WTURN6',wturn6,1.0D0)
118       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
119       call reada(weightcard,'WBOND',wbond,1.0D0)
120       call reada(weightcard,'WTOR',wtor,1.0D0)
121       call reada(weightcard,'WTORD',wtor_d,1.0D0)
122       call reada(weightcard,'WANG',wang,1.0D0)
123       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
124       call reada(weightcard,'SCAL14',scal14,0.4D0)
125       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
126       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
127       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
128       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
129       if (index(weightcard,'SOFT').gt.0) ipot=6
130 C 12/1/95 Added weight for the multi-body term WCORR
131       call reada(weightcard,'WCORRH',wcorr,1.0D0)
132       if (wcorr4.gt.0.0d0) wcorr=wcorr4
133       weights(1)=wsc
134       weights(2)=wscp
135       weights(3)=welec
136       weights(4)=wcorr
137       weights(5)=wcorr5
138       weights(6)=wcorr6
139       weights(7)=wel_loc
140       weights(8)=wturn3
141       weights(9)=wturn4
142       weights(10)=wturn6
143       weights(11)=wang
144       weights(12)=wscloc
145       weights(13)=wtor
146       weights(14)=wtor_d
147       weights(15)=wstrain
148       weights(16)=wvdwpp
149       weights(17)=wbond
150       weights(18)=scal14
151       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
152      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
153      &  wturn4,wturn6,wsccor
154    10 format (/'Energy-term weights (unscaled):'//
155      & 'WSCC=   ',f10.6,' (SC-SC)'/
156      & 'WSCP=   ',f10.6,' (SC-p)'/
157      & 'WELEC=  ',f10.6,' (p-p electr)'/
158      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
159      & 'WBOND=  ',f10.6,' (stretching)'/
160      & 'WANG=   ',f10.6,' (bending)'/
161      & 'WSCLOC= ',f10.6,' (SC local)'/
162      & 'WTOR=   ',f10.6,' (torsional)'/
163      & 'WTORD=  ',f10.6,' (double torsional)'/
164      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
165      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
166      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
167      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
168      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
169      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
170      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
171      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
172      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
173       if (wcorr4.gt.0.0d0) then
174         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
175      &   'between contact pairs of peptide groups'
176         write (iout,'(2(a,f5.3/))')
177      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
178      &  'Range of quenching the correlation terms:',2*delt_corr
179       else if (wcorr.gt.0.0d0) then
180         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
181      &   'between contact pairs of peptide groups'
182       endif
183       write (iout,'(a,f8.3)')
184      &  'Scaling factor of 1,4 SC-p interactions:',scal14
185       write (iout,'(a,f8.3)')
186      &  'General scaling factor of SC-p interactions:',scalscp
187       r0_corr=cutoff_corr-delt_corr
188       do i=1,20
189         aad(i,1)=scalscp*aad(i,1)
190         aad(i,2)=scalscp*aad(i,2)
191         bad(i,1)=scalscp*bad(i,1)
192         bad(i,2)=scalscp*bad(i,2)
193       enddo
194
195       call flush(iout)
196       print *,'indpdb=',indpdb,' pdbref=',pdbref
197
198 C Read sequence if not taken from the pdb file.
199       if (iscode.gt.0) then
200         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
201       else
202         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
203       endif
204 C Convert sequence to numeric code
205       do i=1,nres
206         itype(i)=rescode(i,sequence(i),iscode)
207       enddo
208       print *,nres
209       print '(20i4)',(itype(i),i=1,nres)
210
211       do i=1,nres
212 #ifdef PROCOR
213         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
214 #else
215         if (itype(i).eq.21) then
216 #endif
217           itel(i)=0
218 #ifdef PROCOR
219         else if (itype(i+1).ne.20) then
220 #else
221         else if (itype(i).ne.20) then
222 #endif
223           itel(i)=1
224         else
225           itel(i)=2
226         endif
227       enddo
228       write (iout,*) "ITEL"
229       do i=1,nres-1
230         write (iout,*) i,itype(i),itel(i)
231       enddo
232
233       print *,'Call Read_Bridge.'
234       call read_bridge
235       if (with_dihed_constr) then
236
237       read (inp,*) ndih_constr
238       if (ndih_constr.gt.0) then
239         read (inp,*) ftors
240         write (iout,*) 'FTORS',ftors
241         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
242         write (iout,*)
243      &   'There are',ndih_constr,' constraints on phi angles.'
244         do i=1,ndih_constr
245           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
246         enddo
247         do i=1,ndih_constr
248           phi0(i)=deg2rad*phi0(i)
249           drange(i)=deg2rad*drange(i)
250         enddo
251       endif
252
253       endif
254       nnt=1
255       nct=nres
256       print *,'NNT=',NNT,' NCT=',NCT
257       if (itype(1).eq.21) nnt=2
258       if (itype(nres).eq.21) nct=nct-1
259       if (nstart.lt.nnt) nstart=nnt
260       if (nend.gt.nct .or. nend.eq.0) nend=nct
261       write (iout,*) "nstart",nstart," nend",nend
262       nres0=nres
263 c      if (pdbref) then
264 c        read(inp,'(a)') pdbfile
265 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
266 c        open(ipdbin,file=pdbfile,status='old',err=33)
267 c        goto 34 
268 c  33    write (iout,'(a)') 'Error opening PDB file.'
269 c        stop
270 c  34    continue
271 c        print *,'Begin reading pdb data'
272 c        call readpdb
273 c        print *,'Finished reading pdb data'
274 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
275 c        do i=1,nres
276 c          itype_pdb(i)=itype(i)
277 c        enddo
278 c        close (ipdbin)
279 c        write (iout,'(a,i3)') 'nsup=',nsup
280 c        nstart_seq=nnt
281 c        if (nsup.le.(nct-nnt+1)) then
282 c          do i=0,nct-nnt+1-nsup
283 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
284 c              nstart_seq=nnt+i
285 c              goto 111
286 c            endif
287 c          enddo
288 c          write (iout,'(a)') 
289 c     &            'Error - sequences to be superposed do not match.'
290 c          stop
291 c        else
292 c          do i=0,nsup-(nct-nnt+1)
293 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
294 c     &      then
295 c              nstart_sup=nstart_sup+i
296 c              nsup=nct-nnt+1
297 c              goto 111
298 c            endif
299 c          enddo 
300 c          write (iout,'(a)') 
301 c     &            'Error - sequences to be superposed do not match.'
302 c        endif
303 c  111   continue
304 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
305 c     &                 ' nstart_seq=',nstart_seq
306 c      endif
307       call init_int_table
308       call setup_var
309       write (iout,*) "molread: REFSTR",refstr
310       if (refstr) then
311         if (.not.pdbref) then
312           call read_angles(inp,*38)
313           goto 39
314    38     write (iout,'(a)') 'Error reading reference structure.'
315 #ifdef MPL
316           call mp_stopall(Error_Msg)
317 #else
318           stop 'Error reading reference structure'
319 #endif
320    39     call chainbuild     
321           nstart_sup=nnt
322           nstart_seq=nnt
323           nsup=nct-nnt+1
324           do i=1,2*nres
325             do j=1,3
326               cref(j,i)=c(j,i)
327             enddo
328           enddo
329         endif
330         call contact(.true.,ncont_ref,icont_ref)
331       endif
332 c Read distance restraints
333       if (constr_dist.gt.0) then
334         call read_dist_constr
335         call hpb_partition
336       endif
337       return
338       end
339 c-----------------------------------------------------------------------------
340       logical function seq_comp(itypea,itypeb,length)
341       implicit none
342       integer length,itypea(length),itypeb(length)
343       integer i
344       do i=1,length
345         if (itypea(i).ne.itypeb(i)) then
346           seq_comp=.false.
347           return
348         endif
349       enddo
350       seq_comp=.true.
351       return
352       end
353 c-----------------------------------------------------------------------------
354       subroutine read_bridge
355 C Read information about disulfide bridges.
356       implicit none
357       include 'DIMENSIONS'
358       include 'COMMON.IOUNITS'
359       include 'COMMON.GEO'
360       include 'COMMON.VAR'
361       include 'COMMON.INTERACT'
362       include 'COMMON.LOCAL'
363       include 'COMMON.NAMES'
364       include 'COMMON.CHAIN'
365       include 'COMMON.FFIELD'
366       include 'COMMON.SBRIDGE'
367       include 'COMMON.HEADER'
368       include 'COMMON.CONTROL'
369       include 'COMMON.TIME1'
370 #ifdef MPL
371       include 'COMMON.INFO'
372 #endif
373       integer i,j
374 C Read bridging residues.
375       read (inp,*) ns,(iss(i),i=1,ns)
376       print *,'ns=',ns
377 C Check whether the specified bridging residues are cystines.
378       do i=1,ns
379         if (itype(iss(i)).ne.1) then
380           write (iout,'(2a,i3,a)') 
381      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
382      &   ' can form a disulfide bridge?!!!'
383           write (*,'(2a,i3,a)') 
384      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
385      &   ' can form a disulfide bridge?!!!'
386 #ifdef MPL
387          call mp_stopall(error_msg)
388 #else
389          stop
390 #endif
391         endif
392       enddo
393 C Read preformed bridges.
394       if (ns.gt.0) then
395       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
396       if (nss.gt.0) then
397         nhpb=nss
398 C Check if the residues involved in bridges are in the specified list of
399 C bridging residues.
400         do i=1,nss
401           do j=1,i-1
402             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
403      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
404               write (iout,'(a,i3,a)') 'Disulfide pair',i,
405      &      ' contains residues present in other pairs.'
406               write (*,'(a,i3,a)') 'Disulfide pair',i,
407      &      ' contains residues present in other pairs.'
408 #ifdef MPL
409               call mp_stopall(error_msg)
410 #else
411               stop 
412 #endif
413             endif
414           enddo
415           do j=1,ns
416             if (ihpb(i).eq.iss(j)) goto 10
417           enddo
418           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
419    10     continue
420           do j=1,ns
421             if (jhpb(i).eq.iss(j)) goto 20
422           enddo
423           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
424    20     continue
425 c          dhpb(i)=dbr
426 c          forcon(i)=fbr
427         enddo
428         do i=1,nss
429           ihpb(i)=ihpb(i)+nres
430           jhpb(i)=jhpb(i)+nres
431         enddo
432       endif
433       endif
434       return
435       end
436 c----------------------------------------------------------------------------
437       subroutine read_angles(kanal,*)
438       implicit none
439       include 'DIMENSIONS'
440       include 'COMMON.GEO'
441       include 'COMMON.VAR'
442       include 'COMMON.CHAIN'
443       include 'COMMON.IOUNITS'
444       integer i,kanal
445       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
446       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
447       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
448       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
449       do i=1,nres
450         theta(i)=deg2rad*theta(i)
451         phi(i)=deg2rad*phi(i)
452         alph(i)=deg2rad*alph(i)
453         omeg(i)=deg2rad*omeg(i)
454       enddo
455       return
456    10 return1
457       end
458 c----------------------------------------------------------------------------
459       subroutine reada(rekord,lancuch,wartosc,default)
460       implicit none
461       character*(*) rekord,lancuch
462       double precision wartosc,default
463       integer ilen,iread
464       external ilen
465       iread=index(rekord,lancuch)
466       if (iread.eq.0) then
467         wartosc=default 
468         return
469       endif   
470       iread=iread+ilen(lancuch)+1
471       read (rekord(iread:),*) wartosc
472       return
473       end
474 c----------------------------------------------------------------------------
475       subroutine multreada(rekord,lancuch,tablica,dim,default)
476       implicit none
477       integer dim,i
478       double precision tablica(dim),default
479       character*(*) rekord,lancuch
480       integer ilen,iread
481       external ilen
482       do i=1,dim
483         tablica(i)=default 
484       enddo
485       iread=index(rekord,lancuch)
486       if (iread.eq.0) return
487       iread=iread+ilen(lancuch)+1
488       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
489    10 return
490       end
491 c----------------------------------------------------------------------------
492       subroutine readi(rekord,lancuch,wartosc,default)
493       implicit none
494       character*(*) rekord,lancuch
495       integer wartosc,default
496       integer ilen,iread
497       external ilen
498       iread=index(rekord,lancuch)
499       if (iread.eq.0) then
500         wartosc=default 
501         return
502       endif   
503       iread=iread+ilen(lancuch)+1
504       read (rekord(iread:),*) wartosc
505       return
506       end
507 c----------------------------------------------------------------------------
508       subroutine multreadi(rekord,lancuch,tablica,dim,default)
509       implicit none
510       integer dim,i
511       integer tablica(dim),default
512       character*(*) rekord,lancuch
513       character*80 aux
514       integer ilen,iread
515       external ilen
516       do i=1,dim
517         tablica(i)=default
518       enddo
519       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
520       if (iread.eq.0) return
521       iread=iread+ilen(lancuch)+1
522       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
523    10 return
524       end
525 c----------------------------------------------------------------------------
526       subroutine card_concat(card)
527       include 'DIMENSIONS'
528       include 'COMMON.IOUNITS'
529       character*(*) card
530       character*80 karta,ucase
531       external ilen
532       read (inp,'(a)') karta
533       karta=ucase(karta)
534       card=' '
535       do while (karta(80:80).eq.'&')
536         card=card(:ilen(card)+1)//karta(:79)
537         read (inp,'(a)') karta
538         karta=ucase(karta)
539       enddo
540       card=card(:ilen(card)+1)//karta
541       return
542       end
543 c----------------------------------------------------------------------------
544       subroutine openunits
545       implicit none
546       include 'DIMENSIONS'    
547 #ifdef MPI
548       include "mpif.h"
549       character*3 liczba
550       include "COMMON.MPI"
551 #endif
552       include 'COMMON.IOUNITS'
553       include 'COMMON.CONTROL'
554       integer lenpre,lenpot,ilen
555       external ilen
556       character*16 cformat,cprint
557       character*16 ucase
558       integer lenint,lenout
559       call getenv('INPUT',prefix)
560       call getenv('OUTPUT',prefout)
561       call getenv('INTIN',prefintin)
562       call getenv('COORD',cformat)
563       call getenv('PRINTCOOR',cprint)
564       call getenv('SCRATCHDIR',scratchdir)
565       from_bx=.true.
566       from_cx=.false.
567       if (index(ucase(cformat),'CX').gt.0) then
568         from_cx=.true.
569         from_bx=.false.
570       endif
571       from_cart=.true.
572       lenpre=ilen(prefix)
573       lenout=ilen(prefout)
574       lenint=ilen(prefintin)
575 C Get the names and open the input files
576       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
577 #ifdef MPI
578       write (liczba,'(bz,i3.3)') me
579       outname=prefout(:lenout)//'_clust.out_'//liczba
580 #else
581       outname=prefout(:lenout)//'_clust.out'
582 #endif
583       if (from_bx) then
584         intinname=prefintin(:lenint)//'.bx'
585       else if (from_cx) then
586         intinname=prefintin(:lenint)//'.cx'
587       else
588         intinname=prefintin(:lenint)//'.int'
589       endif
590       rmsname=prefintin(:lenint)//'.rms'
591       open (jplot,file=prefout(:ilen(prefout))//'.tex',
592      &       status='unknown')
593       open (jrms,file=rmsname,status='unknown')
594       open(iout,file=outname,status='unknown')
595 C Get parameter filenames and open the parameter files.
596       call getenv('BONDPAR',bondname)
597       open (ibond,file=bondname,status='old')
598       call getenv('THETPAR',thetname)
599       open (ithep,file=thetname,status='old')
600       call getenv('ROTPAR',rotname)
601       open (irotam,file=rotname,status='old')
602       call getenv('TORPAR',torname)
603       open (itorp,file=torname,status='old')
604       call getenv('TORDPAR',tordname)
605       open (itordp,file=tordname,status='old')
606       call getenv('FOURIER',fouriername)
607       open (ifourier,file=fouriername,status='old')
608       call getenv('ELEPAR',elename)
609       open (ielep,file=elename,status='old')
610       call getenv('SIDEPAR',sidename)
611       open (isidep,file=sidename,status='old')
612       call getenv('SIDEP',sidepname)
613       open (isidep1,file=sidepname,status="old")
614       call getenv('SCCORPAR',sccorname)
615       open (isccor,file=sccorname,status="old")
616 #ifndef OLDSCP
617 C
618 C 8/9/01 In the newest version SCp interaction constants are read from a file
619 C Use -DOLDSCP to use hard-coded constants instead.
620 C
621       call getenv('SCPPAR',scpname)
622       open (iscpp,file=scpname,status='old')
623 #endif
624       return
625       end
626 c-------------------------------------------------------------------------------
627       subroutine read_dist_constr
628       implicit real*8 (a-h,o-z)
629       include 'DIMENSIONS'
630       include 'COMMON.CONTROL'
631       include 'COMMON.CHAIN'
632       include 'COMMON.IOUNITS'
633       include 'COMMON.SBRIDGE'
634       integer ifrag_(2,100),ipair_(2,100)
635       double precision wfrag_(100),wpair_(100)
636       character*500 controlcard
637 c      write (iout,*) "Calling read_dist_constr"
638 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
639       call card_concat(controlcard)
640 c      call flush(iout)
641 c      write (iout,'(a)') controlcard
642       call readi(controlcard,"NFRAG",nfrag_,0)
643       call readi(controlcard,"NPAIR",npair_,0)
644       call readi(controlcard,"NDIST",ndist_,0)
645       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
646       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
647       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
648       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
649       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
650       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
651       write (iout,*) "IFRAG"
652       do i=1,nfrag_
653         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
654       enddo
655       write (iout,*) "IPAIR"
656       do i=1,npair_
657         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
658       enddo
659       call flush(iout)
660       if (.not.refstr .and. nfrag_.gt.0) then
661         write (iout,*) 
662      &  "ERROR: no reference structure to compute distance restraints"
663         write (iout,*)
664      &  "Restraints must be specified explicitly (NDIST=number)"
665         stop 
666       endif
667       if (nfrag_.lt.2 .and. npair_.gt.0) then 
668         write (iout,*) "ERROR: Less than 2 fragments specified",
669      &   " but distance restraints between pairs requested"
670         stop 
671       endif 
672       call flush(iout)
673       do i=1,nfrag_
674         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
675         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
676      &    ifrag_(2,i)=nstart_sup+nsup-1
677 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
678         call flush(iout)
679         if (wfrag_(i).gt.0.0d0) then
680         do j=ifrag_(1,i),ifrag_(2,i)-1
681           do k=j+1,ifrag_(2,i)
682             write (iout,*) "j",j," k",k
683             ddjk=dist(j,k)
684             if (constr_dist.eq.1) then
685             nhpb=nhpb+1
686             ihpb(nhpb)=j
687             jhpb(nhpb)=k
688               dhpb(nhpb)=ddjk
689             forcon(nhpb)=wfrag_(i) 
690             else if (constr_dist.eq.2) then
691               if (ddjk.le.dist_cut) then
692                 nhpb=nhpb+1
693                 ihpb(nhpb)=j
694                 jhpb(nhpb)=k
695                 dhpb(nhpb)=ddjk
696                 forcon(nhpb)=wfrag_(i) 
697               endif
698             else
699               nhpb=nhpb+1
700               ihpb(nhpb)=j
701               jhpb(nhpb)=k
702               dhpb(nhpb)=ddjk
703               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
704             endif
705             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
706      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
707           enddo
708         enddo
709         endif
710       enddo
711       do i=1,npair_
712         if (wpair_(i).gt.0.0d0) then
713         ii = ipair_(1,i)
714         jj = ipair_(2,i)
715         if (ii.gt.jj) then
716           itemp=ii
717           ii=jj
718           jj=itemp
719         endif
720         do j=ifrag_(1,ii),ifrag_(2,ii)
721           do k=ifrag_(1,jj),ifrag_(2,jj)
722             nhpb=nhpb+1
723             ihpb(nhpb)=j
724             jhpb(nhpb)=k
725             forcon(nhpb)=wpair_(i)
726             dhpb(nhpb)=dist(j,k)
727             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
728      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
729           enddo
730         enddo
731         endif
732       enddo 
733       do i=1,ndist_
734         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
735      &     ibecarb(i),forcon(nhpb+1)
736         if (forcon(nhpb+1).gt.0.0d0) then
737           nhpb=nhpb+1
738           if (ibecarb(i).gt.0) then
739             ihpb(i)=ihpb(i)+nres
740             jhpb(i)=jhpb(i)+nres
741           endif
742           if (dhpb(nhpb).eq.0.0d0) 
743      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
744         endif
745       enddo
746       do i=1,nhpb
747           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
748      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
749       enddo
750       call flush(iout)
751       return
752       end