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 character*320 controlcard,ucase
23 read (INP,'(a80)') titel
24 call card_concat(controlcard)
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 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
34 iscode=index(controlcard,'ONE_LETTER')
35 tree=(index(controlcard,'MAKE_TREE').gt.0)
36 min_var=(index(controlcard,'MINVAR').gt.0)
37 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
38 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
39 call readi(controlcard,'NCUT',ncut,0)
40 call readi(controlcard,'NCLUST',nclust,5)
41 call readi(controlcard,'NSTART',nstart,0)
42 call readi(controlcard,'NEND',nend,0)
43 call reada(controlcard,'ECUT',ecut,10.0d0)
44 call reada(controlcard,'PROB',prob_limit,0.99d0)
45 write (iout,*) "Probability limit",prob_limit
46 lgrp=(index(controlcard,'LGRP').gt.0)
47 caonly=(index(controlcard,'CA_ONLY').gt.0)
48 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
50 & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
51 call readi(controlcard,'IOPT',iopt,2)
52 lside = index(controlcard,"SIDE").gt.0
53 efree = index(controlcard,"EFREE").gt.0
54 call readi(controlcard,'NTEMP',nT,1)
55 write (iout,*) "nT",nT
56 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
57 write (iout,*) "nT",nT
58 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
60 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
62 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
63 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
64 lprint_int=index(controlcard,"PRINT_INT") .gt.0
65 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
66 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
67 write (iout,*) "with_dihed_constr ",with_dihed_constr,
68 & " CONSTR_DIST",constr_dist
73 c--------------------------------------------------------------------------
76 C Read molecular data.
80 include 'COMMON.IOUNITS'
83 include 'COMMON.INTERACT'
84 include 'COMMON.LOCAL'
85 include 'COMMON.NAMES'
86 include 'COMMON.CHAIN'
87 include 'COMMON.FFIELD'
88 include 'COMMON.SBRIDGE'
89 include 'COMMON.HEADER'
90 include 'COMMON.CONTROL'
91 include 'COMMON.CONTACTS'
92 include 'COMMON.TIME1'
93 include 'COMMON.TORCNSTR'
97 character*4 sequence(maxres)
98 character*800 weightcard
100 double precision x(maxvar)
101 integer itype_pdb(maxres)
107 C Read weights of the subsequent energy terms.
108 call card_concat(weightcard)
109 call reada(weightcard,'WSC',wsc,1.0d0)
110 call reada(weightcard,'WLONG',wsc,wsc)
111 call reada(weightcard,'WSCP',wscp,1.0d0)
112 call reada(weightcard,'WELEC',welec,1.0D0)
113 call reada(weightcard,'WVDWPP',wvdwpp,welec)
114 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
115 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
116 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
117 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
118 call reada(weightcard,'WTURN3',wturn3,1.0D0)
119 call reada(weightcard,'WTURN4',wturn4,1.0D0)
120 call reada(weightcard,'WTURN6',wturn6,1.0D0)
121 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
122 call reada(weightcard,'WBOND',wbond,1.0D0)
123 call reada(weightcard,'WTOR',wtor,1.0D0)
124 call reada(weightcard,'WTORD',wtor_d,1.0D0)
125 call reada(weightcard,'WANG',wang,1.0D0)
126 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
127 call reada(weightcard,'SCAL14',scal14,0.4D0)
128 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
129 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
130 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
131 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
132 call reada(weightcard,"D0CM",d0cm,3.78d0)
133 call reada(weightcard,"AKCM",akcm,15.1d0)
134 call reada(weightcard,"AKTH",akth,11.0d0)
135 call reada(weightcard,"AKCT",akct,12.0d0)
136 call reada(weightcard,"V1SS",v1ss,-1.08d0)
137 call reada(weightcard,"V2SS",v2ss,7.61d0)
138 call reada(weightcard,"V3SS",v3ss,13.7d0)
139 call reada(weightcard,"EBR",ebr,-5.50D0)
140 if (index(weightcard,'SOFT').gt.0) ipot=6
141 C 12/1/95 Added weight for the multi-body term WCORR
142 call reada(weightcard,'WCORRH',wcorr,1.0D0)
145 dyn_ssbond_ij(i,j)=1.0d300
148 call reada(weightcard,"HT",Ht,0.0D0)
150 ss_depth=ebr/wsc-0.25*eps(1,1)
151 Ht=Ht/wsc-0.25*eps(1,1)
152 akcm=akcm*wstrain/wsc
153 akth=akth*wstrain/wsc
154 akct=akct*wstrain/wsc
155 v1ss=v1ss*wstrain/wsc
156 v2ss=v2ss*wstrain/wsc
157 v3ss=v3ss*wstrain/wsc
159 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
161 write (iout,'(/a)') "Disulfide bridge parameters:"
162 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
163 write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
164 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
165 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
166 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
168 write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
169 if (wcorr4.gt.0.0d0) wcorr=wcorr4
188 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
189 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
190 & wturn4,wturn6,wsccor
191 10 format (/'Energy-term weights (unscaled):'//
192 & 'WSCC= ',f10.6,' (SC-SC)'/
193 & 'WSCP= ',f10.6,' (SC-p)'/
194 & 'WELEC= ',f10.6,' (p-p electr)'/
195 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
196 & 'WBOND= ',f10.6,' (stretching)'/
197 & 'WANG= ',f10.6,' (bending)'/
198 & 'WSCLOC= ',f10.6,' (SC local)'/
199 & 'WTOR= ',f10.6,' (torsional)'/
200 & 'WTORD= ',f10.6,' (double torsional)'/
201 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
202 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
203 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
204 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
205 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
206 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
207 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
208 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
209 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
210 if (wcorr4.gt.0.0d0) then
211 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
212 & 'between contact pairs of peptide groups'
213 write (iout,'(2(a,f5.3/))')
214 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
215 & 'Range of quenching the correlation terms:',2*delt_corr
216 else if (wcorr.gt.0.0d0) then
217 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
218 & 'between contact pairs of peptide groups'
220 write (iout,'(a,f8.3)')
221 & 'Scaling factor of 1,4 SC-p interactions:',scal14
222 write (iout,'(a,f8.3)')
223 & 'General scaling factor of SC-p interactions:',scalscp
224 r0_corr=cutoff_corr-delt_corr
226 aad(i,1)=scalscp*aad(i,1)
227 aad(i,2)=scalscp*aad(i,2)
228 bad(i,1)=scalscp*bad(i,1)
229 bad(i,2)=scalscp*bad(i,2)
233 print *,'indpdb=',indpdb,' pdbref=',pdbref
235 C Read sequence if not taken from the pdb file.
236 if (iscode.gt.0) then
237 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
239 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
241 C Convert sequence to numeric code
243 itype(i)=rescode(i,sequence(i),iscode)
245 if (itype(2).eq.10) then
247 & "Glycine is the first full residue, initial dummy deleted"
253 if (itype(nres).eq.10) then
255 & "Glycine is the last full residue, terminal dummy deleted"
259 print '(20i4)',(itype(i),i=1,nres)
263 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
265 if (itype(i).eq.21) then
269 else if (itype(i+1).ne.20) then
271 else if (itype(i).ne.20) then
278 write (iout,*) "ITEL"
280 write (iout,*) i,itype(i),itel(i)
283 print *,'Call Read_Bridge.'
285 if (with_dihed_constr) then
287 read (inp,*) ndih_constr
288 if (ndih_constr.gt.0) then
290 write (iout,*) 'FTORS',ftors
291 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
293 & 'There are',ndih_constr,' constraints on phi angles.'
295 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
298 phi0(i)=deg2rad*phi0(i)
299 drange(i)=deg2rad*drange(i)
306 print *,'NNT=',NNT,' NCT=',NCT
307 if (itype(1).eq.21) nnt=2
308 if (itype(nres).eq.21) nct=nct-1
309 if (nstart.lt.nnt) nstart=nnt
310 if (nend.gt.nct .or. nend.eq.0) nend=nct
311 write (iout,*) "nstart",nstart," nend",nend
314 c read(inp,'(a)') pdbfile
315 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
316 c open(ipdbin,file=pdbfile,status='old',err=33)
318 c 33 write (iout,'(a)') 'Error opening PDB file.'
321 c print *,'Begin reading pdb data'
323 c print *,'Finished reading pdb data'
324 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
326 c itype_pdb(i)=itype(i)
329 c write (iout,'(a,i3)') 'nsup=',nsup
331 c if (nsup.le.(nct-nnt+1)) then
332 c do i=0,nct-nnt+1-nsup
333 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
339 c & 'Error - sequences to be superposed do not match.'
342 c do i=0,nsup-(nct-nnt+1)
343 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
345 c nstart_sup=nstart_sup+i
351 c & 'Error - sequences to be superposed do not match.'
354 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
355 c & ' nstart_seq=',nstart_seq
359 write (iout,*) "molread: REFSTR",refstr
361 if (.not.pdbref) then
362 call read_angles(inp,*38)
364 38 write (iout,'(a)') 'Error reading reference structure.'
366 call mp_stopall(Error_Msg)
368 stop 'Error reading reference structure'
380 call contact(print_contact_map,ncont_ref,icont_ref)
382 c Read distance restraints
383 if (constr_dist.gt.0) then
384 call read_dist_constr
389 c-----------------------------------------------------------------------------
390 logical function seq_comp(itypea,itypeb,length)
392 integer length,itypea(length),itypeb(length)
395 if (itypea(i).ne.itypeb(i)) then
403 c-----------------------------------------------------------------------------
404 subroutine read_bridge
405 C Read information about disulfide bridges.
408 include 'COMMON.IOUNITS'
411 include 'COMMON.INTERACT'
412 include 'COMMON.LOCAL'
413 include 'COMMON.NAMES'
414 include 'COMMON.CHAIN'
415 include 'COMMON.FFIELD'
416 include 'COMMON.SBRIDGE'
417 include 'COMMON.HEADER'
418 include 'COMMON.CONTROL'
419 include 'COMMON.TIME1'
421 include 'COMMON.INFO'
424 C Read bridging residues.
425 read (inp,*) ns,(iss(i),i=1,ns)
427 C Check whether the specified bridging residues are cystines.
429 if (itype(iss(i)).ne.1) then
430 write (iout,'(2a,i3,a)')
431 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
432 & ' can form a disulfide bridge?!!!'
433 write (*,'(2a,i3,a)')
434 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
435 & ' can form a disulfide bridge?!!!'
437 call mp_stopall(error_msg)
443 C Read preformed bridges.
445 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
448 C Check if the residues involved in bridges are in the specified list of
452 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
453 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
454 write (iout,'(a,i3,a)') 'Disulfide pair',i,
455 & ' contains residues present in other pairs.'
456 write (*,'(a,i3,a)') 'Disulfide pair',i,
457 & ' contains residues present in other pairs.'
459 call mp_stopall(error_msg)
466 if (ihpb(i).eq.iss(j)) goto 10
468 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
471 if (jhpb(i).eq.iss(j)) goto 20
473 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
484 if (ns.gt.0.and.dyn_ss) then
488 forcon(i-nss)=forcon(i)
495 dyn_ss_mask(iss(i))=.true.
496 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
499 print *, "Leaving brigde read"
502 c----------------------------------------------------------------------------
503 subroutine read_angles(kanal,*)
508 include 'COMMON.CHAIN'
509 include 'COMMON.IOUNITS'
511 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
512 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
513 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
514 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
516 theta(i)=deg2rad*theta(i)
517 phi(i)=deg2rad*phi(i)
518 alph(i)=deg2rad*alph(i)
519 omeg(i)=deg2rad*omeg(i)
524 c----------------------------------------------------------------------------
525 subroutine reada(rekord,lancuch,wartosc,default)
527 character*(*) rekord,lancuch
528 double precision wartosc,default
531 iread=index(rekord,lancuch)
536 iread=iread+ilen(lancuch)+1
537 read (rekord(iread:),*) wartosc
540 c----------------------------------------------------------------------------
541 subroutine multreada(rekord,lancuch,tablica,dim,default)
544 double precision tablica(dim),default
545 character*(*) rekord,lancuch
551 iread=index(rekord,lancuch)
552 if (iread.eq.0) return
553 iread=iread+ilen(lancuch)+1
554 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
557 c----------------------------------------------------------------------------
558 subroutine readi(rekord,lancuch,wartosc,default)
560 character*(*) rekord,lancuch
561 integer wartosc,default
564 iread=index(rekord,lancuch)
569 iread=iread+ilen(lancuch)+1
570 read (rekord(iread:),*) wartosc
573 c----------------------------------------------------------------------------
574 subroutine multreadi(rekord,lancuch,tablica,dim,default)
577 integer tablica(dim),default
578 character*(*) rekord,lancuch
585 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
586 if (iread.eq.0) return
587 iread=iread+ilen(lancuch)+1
588 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
591 c----------------------------------------------------------------------------
592 subroutine card_concat(card)
594 include 'COMMON.IOUNITS'
596 character*80 karta,ucase
598 read (inp,'(a)') karta
601 do while (karta(80:80).eq.'&')
602 card=card(:ilen(card)+1)//karta(:79)
603 read (inp,'(a)') karta
606 card=card(:ilen(card)+1)//karta
609 c----------------------------------------------------------------------------
618 include 'COMMON.IOUNITS'
619 include 'COMMON.CONTROL'
620 integer lenpre,lenpot,ilen
622 character*16 cformat,cprint
624 integer lenint,lenout
625 call getenv('INPUT',prefix)
626 call getenv('OUTPUT',prefout)
627 call getenv('INTIN',prefintin)
628 call getenv('COORD',cformat)
629 call getenv('PRINTCOOR',cprint)
630 call getenv('SCRATCHDIR',scratchdir)
633 if (index(ucase(cformat),'CX').gt.0) then
640 lenint=ilen(prefintin)
641 C Get the names and open the input files
642 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
644 write (liczba,'(bz,i3.3)') me
645 outname=prefout(:lenout)//'_clust.out_'//liczba
647 outname=prefout(:lenout)//'_clust.out'
650 intinname=prefintin(:lenint)//'.bx'
651 else if (from_cx) then
652 intinname=prefintin(:lenint)//'.cx'
654 intinname=prefintin(:lenint)//'.int'
656 rmsname=prefintin(:lenint)//'.rms'
657 open (jplot,file=prefout(:ilen(prefout))//'.tex',
659 open (jrms,file=rmsname,status='unknown')
660 open(iout,file=outname,status='unknown')
661 C Get parameter filenames and open the parameter files.
662 call getenv('BONDPAR',bondname)
663 open (ibond,file=bondname,status='old')
664 call getenv('THETPAR',thetname)
665 open (ithep,file=thetname,status='old')
666 call getenv('ROTPAR',rotname)
667 open (irotam,file=rotname,status='old')
668 call getenv('TORPAR',torname)
669 open (itorp,file=torname,status='old')
670 call getenv('TORDPAR',tordname)
671 open (itordp,file=tordname,status='old')
672 call getenv('FOURIER',fouriername)
673 open (ifourier,file=fouriername,status='old')
674 call getenv('ELEPAR',elename)
675 open (ielep,file=elename,status='old')
676 call getenv('SIDEPAR',sidename)
677 open (isidep,file=sidename,status='old')
678 call getenv('SIDEP',sidepname)
679 open (isidep1,file=sidepname,status="old")
680 call getenv('SCCORPAR',sccorname)
681 open (isccor,file=sccorname,status="old")
684 C 8/9/01 In the newest version SCp interaction constants are read from a file
685 C Use -DOLDSCP to use hard-coded constants instead.
687 call getenv('SCPPAR',scpname)
688 open (iscpp,file=scpname,status='old')
692 c-------------------------------------------------------------------------------
693 subroutine read_dist_constr
694 implicit real*8 (a-h,o-z)
696 include 'COMMON.CONTROL'
697 include 'COMMON.CHAIN'
698 include 'COMMON.IOUNITS'
699 include 'COMMON.SBRIDGE'
700 integer ifrag_(2,100),ipair_(2,100)
701 double precision wfrag_(100),wpair_(100)
702 character*500 controlcard
703 c write (iout,*) "Calling read_dist_constr"
704 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
705 call card_concat(controlcard)
707 c write (iout,'(a)') controlcard
708 call readi(controlcard,"NFRAG",nfrag_,0)
709 call readi(controlcard,"NPAIR",npair_,0)
710 call readi(controlcard,"NDIST",ndist_,0)
711 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
712 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
713 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
714 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
715 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
716 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
717 write (iout,*) "IFRAG"
719 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
721 write (iout,*) "IPAIR"
723 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
726 if (.not.refstr .and. nfrag_.gt.0) then
728 & "ERROR: no reference structure to compute distance restraints"
730 & "Restraints must be specified explicitly (NDIST=number)"
733 if (nfrag_.lt.2 .and. npair_.gt.0) then
734 write (iout,*) "ERROR: Less than 2 fragments specified",
735 & " but distance restraints between pairs requested"
740 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
741 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
742 & ifrag_(2,i)=nstart_sup+nsup-1
743 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
745 if (wfrag_(i).gt.0.0d0) then
746 do j=ifrag_(1,i),ifrag_(2,i)-1
748 write (iout,*) "j",j," k",k
750 if (constr_dist.eq.1) then
755 forcon(nhpb)=wfrag_(i)
756 else if (constr_dist.eq.2) then
757 if (ddjk.le.dist_cut) then
762 forcon(nhpb)=wfrag_(i)
769 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
771 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
772 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
778 if (wpair_(i).gt.0.0d0) then
786 do j=ifrag_(1,ii),ifrag_(2,ii)
787 do k=ifrag_(1,jj),ifrag_(2,jj)
791 forcon(nhpb)=wpair_(i)
793 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
794 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
800 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
801 & ibecarb(i),forcon(nhpb+1)
802 if (forcon(nhpb+1).gt.0.0d0) then
804 if (ibecarb(i).gt.0) then
808 if (dhpb(nhpb).eq.0.0d0)
809 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
813 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
814 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)