1 subroutine read_control
8 include 'COMMON.IOUNITS'
10 include 'COMMON.SBRIDGE'
11 include 'COMMON.CONTROL'
12 include 'COMMON.CLUSTER'
13 include 'COMMON.CHAIN'
14 include 'COMMON.HEADER'
15 include 'COMMON.FFIELD'
17 include 'COMMON.INTERACT'
18 include "COMMON.SPLITELE"
19 character*320 controlcard,ucase
23 integer i,i1,i2,it1,it2
25 read (INP,'(a80)') titel
26 call card_concat(controlcard)
28 call readi(controlcard,'NRES',nres,0)
29 call readi(controlcard,'RESCALE',rescale_mode,2)
30 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
31 write (iout,*) "DISTCHAINMAX",distchainmax
32 C Reading the dimensions of box in x,y,z coordinates
33 call reada(controlcard,'BOXX',boxxsize,100.0d0)
34 call reada(controlcard,'BOXY',boxysize,100.0d0)
35 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
36 c Cutoff range for interactions
37 call reada(controlcard,"R_CUT",r_cut,15.0d0)
38 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
39 call readi(controlcard,'PDBOUT',outpdb,0)
40 call readi(controlcard,'MOL2OUT',outmol2,0)
41 refstr=(index(controlcard,'REFSTR').gt.0)
42 write (iout,*) "REFSTR",refstr
43 pdbref=(index(controlcard,'PDBREF').gt.0)
44 iscode=index(controlcard,'ONE_LETTER')
45 tree=(index(controlcard,'MAKE_TREE').gt.0)
46 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
47 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
48 write (iout,*) "with_dihed_constr ",with_dihed_constr,
49 & " CONSTR_DIST",constr_dist
50 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
51 write (iout,*) "with_theta_constr ",with_theta_constr
53 min_var=(index(controlcard,'MINVAR').gt.0)
54 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
55 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
56 call readi(controlcard,'NCUT',ncut,1)
57 call readi(controlcard,'SYM',symetr,1)
58 write (iout,*) 'sym', symetr
59 call readi(controlcard,'NSTART',nstart,0)
60 call readi(controlcard,'NEND',nend,0)
61 call reada(controlcard,'ECUT',ecut,10.0d0)
62 call reada(controlcard,'PROB',prob_limit,0.99d0)
63 write (iout,*) "Probability limit",prob_limit
64 lgrp=(index(controlcard,'LGRP').gt.0)
65 caonly=(index(controlcard,'CA_ONLY').gt.0)
66 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
67 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
68 call readi(controlcard,'IOPT',iopt,2)
69 lside = index(controlcard,"SIDE").gt.0
70 efree = index(controlcard,"EFREE").gt.0
71 call readi(controlcard,'NTEMP',nT,1)
72 write (iout,*) "nT",nT
73 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
74 write (iout,*) "nT",nT
75 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
77 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
79 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
80 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
81 lprint_int=index(controlcard,"PRINT_INT") .gt.0
85 c--------------------------------------------------------------------------
88 C Read molecular data.
92 include 'COMMON.IOUNITS'
95 include 'COMMON.INTERACT'
96 include 'COMMON.LOCAL'
97 include 'COMMON.NAMES'
98 include 'COMMON.CHAIN'
99 include 'COMMON.FFIELD'
100 include 'COMMON.SBRIDGE'
101 include 'COMMON.HEADER'
102 include 'COMMON.CONTROL'
103 include 'COMMON.CONTACTS'
104 include 'COMMON.TIME1'
105 include 'COMMON.TORCNSTR'
107 include 'COMMON.INFO'
109 character*4 sequence(maxres)
110 character*800 weightcard
112 double precision x(maxvar)
113 integer itype_pdb(maxres)
115 integer i,j,kkk,i1,i2,it1,it2
119 C Read weights of the subsequent energy terms.
120 call card_concat(weightcard)
121 call reada(weightcard,'WSC',wsc,1.0d0)
122 call reada(weightcard,'WLONG',wsc,wsc)
123 call reada(weightcard,'WSCP',wscp,1.0d0)
124 call reada(weightcard,'WELEC',welec,1.0D0)
125 call reada(weightcard,'WVDWPP',wvdwpp,welec)
126 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
127 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
128 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
129 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
130 call reada(weightcard,'WTURN3',wturn3,1.0D0)
131 call reada(weightcard,'WTURN4',wturn4,1.0D0)
132 call reada(weightcard,'WTURN6',wturn6,1.0D0)
133 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
134 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
135 call reada(weightcard,'WBOND',wbond,1.0D0)
136 call reada(weightcard,'WTOR',wtor,1.0D0)
137 call reada(weightcard,'WTORD',wtor_d,1.0D0)
138 call reada(weightcard,'WANG',wang,1.0D0)
139 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
140 call reada(weightcard,'SCAL14',scal14,0.4D0)
141 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
142 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
143 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
144 if (index(weightcard,'SOFT').gt.0) ipot=6
145 call reada(weightcard,"D0CM",d0cm,3.78d0)
146 call reada(weightcard,"AKCM",akcm,15.1d0)
147 call reada(weightcard,"AKTH",akth,11.0d0)
148 call reada(weightcard,"AKCT",akct,12.0d0)
149 call reada(weightcard,"V1SS",v1ss,-1.08d0)
150 call reada(weightcard,"V2SS",v2ss,7.61d0)
151 call reada(weightcard,"V3SS",v3ss,13.7d0)
152 call reada(weightcard,"EBR",ebr,-5.50D0)
153 call reada(weightcard,"ATRISS",atriss,0.301D0)
154 call reada(weightcard,"BTRISS",btriss,0.021D0)
155 call reada(weightcard,"CTRISS",ctriss,1.001D0)
156 call reada(weightcard,"DTRISS",dtriss,1.001D0)
157 write (iout,*) "ATRISS=", atriss
158 write (iout,*) "BTRISS=", btriss
159 write (iout,*) "CTRISS=", ctriss
160 write (iout,*) "DTRISS=", dtriss
161 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
163 dyn_ss_mask(i)=.false.
167 dyn_ssbond_ij(i,j)=1.0d300
170 call reada(weightcard,"HT",Ht,0.0D0)
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
181 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
183 write (iout,'(/a)') "Disulfide bridge parameters:"
184 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
185 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
186 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
187 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
190 C 12/1/95 Added weight for the multi-body term WCORR
191 call reada(weightcard,'WCORRH',wcorr,1.0D0)
192 if (wcorr4.gt.0.0d0) wcorr=wcorr4
212 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
213 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
214 & wturn4,wturn6,wsccor
215 10 format (/'Energy-term weights (unscaled):'//
216 & 'WSCC= ',f10.6,' (SC-SC)'/
217 & 'WSCP= ',f10.6,' (SC-p)'/
218 & 'WELEC= ',f10.6,' (p-p electr)'/
219 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
220 & 'WBOND= ',f10.6,' (stretching)'/
221 & 'WANG= ',f10.6,' (bending)'/
222 & 'WSCLOC= ',f10.6,' (SC local)'/
223 & 'WTOR= ',f10.6,' (torsional)'/
224 & 'WTORD= ',f10.6,' (double torsional)'/
225 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
226 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
227 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
228 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
229 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
230 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
231 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
232 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
233 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
235 if (wcorr4.gt.0.0d0) then
236 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
237 & 'between contact pairs of peptide groups'
238 write (iout,'(2(a,f5.3/))')
239 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
240 & 'Range of quenching the correlation terms:',2*delt_corr
241 else if (wcorr.gt.0.0d0) then
242 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
243 & 'between contact pairs of peptide groups'
245 write (iout,'(a,f8.3)')
246 & 'Scaling factor of 1,4 SC-p interactions:',scal14
247 write (iout,'(a,f8.3)')
248 & 'General scaling factor of SC-p interactions:',scalscp
249 r0_corr=cutoff_corr-delt_corr
251 aad(i,1)=scalscp*aad(i,1)
252 aad(i,2)=scalscp*aad(i,2)
253 bad(i,1)=scalscp*bad(i,1)
254 bad(i,2)=scalscp*bad(i,2)
258 print *,'indpdb=',indpdb,' pdbref=',pdbref
260 C Read sequence if not taken from the pdb file.
261 if (iscode.gt.0) then
262 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
264 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
266 C Convert sequence to numeric code
268 itype(i)=rescode(i,sequence(i),iscode)
271 print '(20i4)',(itype(i),i=1,nres)
275 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
277 if (itype(i).eq.ntyp1) then
281 else if (iabs(itype(i+1)).ne.20) then
283 else if (iabs(itype(i)).ne.20) then
290 write (iout,*) "ITEL"
292 write (iout,*) i,itype(i),itel(i)
295 print *,'Call Read_Bridge.'
297 C this fragment reads diheadral constrains
298 if (with_dihed_constr) then
300 read (inp,*) ndih_constr
301 if (ndih_constr.gt.0) then
303 C write (iout,*) 'FTORS',ftors
304 C ftors is the force constant for torsional quartic constrains
305 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
308 & 'There are',ndih_constr,' constraints on phi angles.'
310 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
314 phi0(i)=deg2rad*phi0(i)
315 drange(i)=deg2rad*drange(i)
317 endif ! endif ndif_constr.gt.0
318 endif ! with_dihed_constr
319 if (with_theta_constr) then
320 C with_theta_constr is keyword allowing for occurance of theta constrains
321 read (inp,*) ntheta_constr
322 C ntheta_constr is the number of theta constrains
323 if (ntheta_constr.gt.0) then
325 read (inp,*) (itheta_constr(i),theta_constr0(i),
326 & theta_drange(i),for_thet_constr(i),
328 C the above code reads from 1 to ntheta_constr
329 C itheta_constr(i) residue i for which is theta_constr
330 C theta_constr0 the global minimum value
331 C theta_drange is range for which there is no energy penalty
332 C for_thet_constr is the force constant for quartic energy penalty
334 C if(me.eq.king.or..not.out1file)then
336 & 'There are',ntheta_constr,' constraints on phi angles.'
338 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
344 theta_constr0(i)=deg2rad*theta_constr0(i)
345 theta_drange(i)=deg2rad*theta_drange(i)
347 C if(me.eq.king.or..not.out1file)
348 C & write (iout,*) 'FTORS',ftors
349 C do i=1,ntheta_constr
350 C ii = itheta_constr(i)
351 C thetabound(1,ii) = phi0(i)-drange(i)
352 C thetabound(2,ii) = phi0(i)+drange(i)
354 endif ! ntheta_constr.gt.0
355 endif! with_theta_constr
359 print *,'NNT=',NNT,' NCT=',NCT
360 if (itype(1).eq.ntyp1) nnt=2
361 if (itype(nres).eq.ntyp1) nct=nct-1
362 if (nstart.lt.nnt) nstart=nnt
363 if (nend.gt.nct .or. nend.eq.0) nend=nct
364 write (iout,*) "nstart",nstart," nend",nend
367 c read(inp,'(a)') pdbfile
368 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
369 c open(ipdbin,file=pdbfile,status='old',err=33)
371 c 33 write (iout,'(a)') 'Error opening PDB file.'
374 c print *,'Begin reading pdb data'
376 c print *,'Finished reading pdb data'
377 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
379 c itype_pdb(i)=itype(i)
382 c write (iout,'(a,i3)') 'nsup=',nsup
384 c if (nsup.le.(nct-nnt+1)) then
385 c do i=0,nct-nnt+1-nsup
386 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
392 c & 'Error - sequences to be superposed do not match.'
395 c do i=0,nsup-(nct-nnt+1)
396 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
398 c nstart_sup=nstart_sup+i
404 c & 'Error - sequences to be superposed do not match.'
407 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
408 c & ' nstart_seq=',nstart_seq
412 write (iout,*) "molread: REFSTR",refstr
414 if (.not.pdbref) then
415 call read_angles(inp,*38)
417 38 write (iout,'(a)') 'Error reading reference structure.'
419 call mp_stopall(Error_Msg)
421 stop 'Error reading reference structure'
434 call contact(.true.,ncont_ref,icont_ref)
437 C write (iout,'(/a,i3,a)')
438 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
439 write (iout,'(20i4)') (iss(i),i=1,ns)
441 write(iout,*)"Running with dynamic disulfide-bond formation"
443 write (iout,'(/a/)') 'Pre-formed links are:'
449 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
450 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
456 if (ns.gt.0.and.dyn_ss) then
460 forcon(i-nss)=forcon(i)
467 dyn_ss_mask(iss(i))=.true.
470 c Read distance restraints
471 if (constr_dist.gt.0) then
472 call read_dist_constr
477 c-----------------------------------------------------------------------------
478 logical function seq_comp(itypea,itypeb,length)
480 integer length,itypea(length),itypeb(length)
483 if (itypea(i).ne.itypeb(i)) then
491 c-----------------------------------------------------------------------------
492 subroutine read_bridge
493 C Read information about disulfide bridges.
496 include 'COMMON.IOUNITS'
499 include 'COMMON.INTERACT'
500 include 'COMMON.LOCAL'
501 include 'COMMON.NAMES'
502 include 'COMMON.CHAIN'
503 include 'COMMON.FFIELD'
504 include 'COMMON.SBRIDGE'
505 include 'COMMON.HEADER'
506 include 'COMMON.CONTROL'
507 include 'COMMON.TIME1'
509 include 'COMMON.INFO'
512 C Read bridging residues.
513 read (inp,*) ns,(iss(i),i=1,ns)
515 C Check whether the specified bridging residues are cystines.
517 if (itype(iss(i)).ne.1) then
518 write (iout,'(2a,i3,a)')
519 & 'Do you REALLY think that the residue ',
520 & restyp(itype(iss(i))),i,
521 & ' can form a disulfide bridge?!!!'
522 write (*,'(2a,i3,a)')
523 & 'Do you REALLY think that the residue ',
524 & restyp(itype(iss(i))),i,
525 & ' can form a disulfide bridge?!!!'
527 call mp_stopall(error_msg)
533 C Read preformed bridges.
535 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
538 C Check if the residues involved in bridges are in the specified list of
542 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
543 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
544 write (iout,'(a,i3,a)') 'Disulfide pair',i,
545 & ' contains residues present in other pairs.'
546 write (*,'(a,i3,a)') 'Disulfide pair',i,
547 & ' contains residues present in other pairs.'
549 call mp_stopall(error_msg)
556 if (ihpb(i).eq.iss(j)) goto 10
558 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
561 if (jhpb(i).eq.iss(j)) goto 20
563 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
576 c----------------------------------------------------------------------------
577 subroutine read_angles(kanal,*)
582 include 'COMMON.CHAIN'
583 include 'COMMON.IOUNITS'
585 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
586 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
587 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
588 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
590 theta(i)=deg2rad*theta(i)
591 phi(i)=deg2rad*phi(i)
592 alph(i)=deg2rad*alph(i)
593 omeg(i)=deg2rad*omeg(i)
598 c----------------------------------------------------------------------------
599 subroutine reada(rekord,lancuch,wartosc,default)
601 character*(*) rekord,lancuch
602 double precision wartosc,default
605 iread=index(rekord,lancuch)
610 iread=iread+ilen(lancuch)+1
611 read (rekord(iread:),*) wartosc
614 c----------------------------------------------------------------------------
615 subroutine multreada(rekord,lancuch,tablica,dim,default)
618 double precision tablica(dim),default
619 character*(*) rekord,lancuch
625 iread=index(rekord,lancuch)
626 if (iread.eq.0) return
627 iread=iread+ilen(lancuch)+1
628 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
631 c----------------------------------------------------------------------------
632 subroutine readi(rekord,lancuch,wartosc,default)
634 character*(*) rekord,lancuch
635 integer wartosc,default
638 iread=index(rekord,lancuch)
643 iread=iread+ilen(lancuch)+1
644 read (rekord(iread:),*) wartosc
647 C----------------------------------------------------------------------
648 subroutine multreadi(rekord,lancuch,tablica,dim,default)
651 integer tablica(dim),default
652 character*(*) rekord,lancuch
659 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
660 if (iread.eq.0) return
661 iread=iread+ilen(lancuch)+1
662 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
666 c----------------------------------------------------------------------------
667 subroutine card_concat(card)
669 include 'COMMON.IOUNITS'
671 character*80 karta,ucase
673 read (inp,'(a)') karta
676 do while (karta(80:80).eq.'&')
677 card=card(:ilen(card)+1)//karta(:79)
678 read (inp,'(a)') karta
681 card=card(:ilen(card)+1)//karta
684 c----------------------------------------------------------------------------
693 include 'COMMON.IOUNITS'
694 include 'COMMON.CONTROL'
695 integer lenpre,lenpot,ilen
697 character*16 cformat,cprint
699 integer lenint,lenout
700 call getenv('INPUT',prefix)
701 call getenv('OUTPUT',prefout)
702 call getenv('INTIN',prefintin)
703 call getenv('COORD',cformat)
704 call getenv('PRINTCOOR',cprint)
705 call getenv('SCRATCHDIR',scratchdir)
708 if (index(ucase(cformat),'CX').gt.0) then
715 lenint=ilen(prefintin)
716 C Get the names and open the input files
717 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
719 write (liczba,'(bz,i3.3)') me
720 outname=prefout(:lenout)//'_clust.out_'//liczba
722 outname=prefout(:lenout)//'_clust.out'
725 intinname=prefintin(:lenint)//'.bx'
726 else if (from_cx) then
727 intinname=prefintin(:lenint)//'.cx'
729 intinname=prefintin(:lenint)//'.int'
731 rmsname=prefintin(:lenint)//'.rms'
732 open (jplot,file=prefout(:ilen(prefout))//'.tex',
734 open (jrms,file=rmsname,status='unknown')
735 open(iout,file=outname,status='unknown')
736 C Get parameter filenames and open the parameter files.
737 call getenv('BONDPAR',bondname)
738 open (ibond,file=bondname,status='old')
739 call getenv('THETPAR',thetname)
740 open (ithep,file=thetname,status='old')
741 call getenv('ROTPAR',rotname)
742 open (irotam,file=rotname,status='old')
743 call getenv('TORPAR',torname)
744 open (itorp,file=torname,status='old')
745 call getenv('TORDPAR',tordname)
746 open (itordp,file=tordname,status='old')
747 call getenv('FOURIER',fouriername)
748 open (ifourier,file=fouriername,status='old')
749 call getenv('ELEPAR',elename)
750 open (ielep,file=elename,status='old')
751 call getenv('SIDEPAR',sidename)
752 open (isidep,file=sidename,status='old')
753 call getenv('SIDEP',sidepname)
754 open (isidep1,file=sidepname,status="old")
755 call getenv('SCCORPAR',sccorname)
756 open (isccor,file=sccorname,status="old")
759 C 8/9/01 In the newest version SCp interaction constants are read from a file
760 C Use -DOLDSCP to use hard-coded constants instead.
762 call getenv('SCPPAR',scpname)
763 open (iscpp,file=scpname,status='old')
767 subroutine read_dist_constr
768 implicit real*8 (a-h,o-z)
773 include 'COMMON.CONTROL'
774 include 'COMMON.CHAIN'
775 include 'COMMON.IOUNITS'
776 include 'COMMON.SBRIDGE'
777 integer ifrag_(2,100),ipair_(2,100)
778 double precision wfrag_(100),wpair_(100)
779 character*500 controlcard
780 logical lprn /.true./
781 write (iout,*) "Calling read_dist_constr"
782 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
784 write(iout,*) "TU sie wywalam?"
785 call card_concat(controlcard)
786 write (iout,*) controlcard
788 call readi(controlcard,"NFRAG",nfrag_,0)
789 call readi(controlcard,"NPAIR",npair_,0)
790 call readi(controlcard,"NDIST",ndist_,0)
791 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
792 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
793 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
794 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
795 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
796 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
797 write (iout,*) "IFRAG"
799 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
801 write (iout,*) "IPAIR"
803 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
807 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
808 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
809 & ifrag_(2,i)=nstart_sup+nsup-1
810 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
812 if (wfrag_(i).gt.0.0d0) then
813 do j=ifrag_(1,i),ifrag_(2,i)-1
815 write (iout,*) "j",j," k",k
817 if (constr_dist.eq.1) then
822 forcon(nhpb)=wfrag_(i)
823 else if (constr_dist.eq.2) then
824 if (ddjk.le.dist_cut) then
829 forcon(nhpb)=wfrag_(i)
836 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
839 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
840 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
846 if (wpair_(i).gt.0.0d0) then
854 do j=ifrag_(1,ii),ifrag_(2,ii)
855 do k=ifrag_(1,jj),ifrag_(2,jj)
859 forcon(nhpb)=wpair_(i)
861 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
862 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
868 if (constr_dist.eq.11) then
869 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
870 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
871 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
872 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
873 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
875 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
877 if (forcon(nhpb+1).gt.0.0d0) then
879 if (ibecarb(i).gt.0) then
883 if (dhpb(nhpb).eq.0.0d0)
884 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
885 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
886 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
887 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)