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,1)
40 call readi(controlcard,'NSTART',nstart,0)
41 call readi(controlcard,'NEND',nend,0)
42 call reada(controlcard,'ECUT',ecut,10.0d0)
43 call reada(controlcard,'PROB',prob_limit,0.99d0)
44 write (iout,*) "Probability limit",prob_limit
45 lgrp=(index(controlcard,'LGRP').gt.0)
46 caonly=(index(controlcard,'CA_ONLY').gt.0)
47 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
48 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
49 call readi(controlcard,'IOPT',iopt,2)
50 lside = index(controlcard,"SIDE").gt.0
51 efree = index(controlcard,"EFREE").gt.0
52 call readi(controlcard,'NTEMP',nT,1)
53 write (iout,*) "nT",nT
54 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
55 write (iout,*) "nT",nT
56 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
58 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
60 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
61 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
62 lprint_int=index(controlcard,"PRINT_INT") .gt.0
63 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
64 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
65 write (iout,*) "with_dihed_constr ",with_dihed_constr,
66 & " CONSTR_DIST",constr_dist
71 c--------------------------------------------------------------------------
74 C Read molecular data.
78 include 'COMMON.IOUNITS'
81 include 'COMMON.INTERACT'
82 include 'COMMON.LOCAL'
83 include 'COMMON.NAMES'
84 include 'COMMON.CHAIN'
85 include 'COMMON.FFIELD'
86 include 'COMMON.SBRIDGE'
87 include 'COMMON.HEADER'
88 include 'COMMON.CONTROL'
89 include 'COMMON.CONTACTS'
90 include 'COMMON.TIME1'
91 include 'COMMON.TORCNSTR'
95 character*4 sequence(maxres)
96 character*800 weightcard
98 double precision x(maxvar)
99 integer itype_pdb(maxres)
105 C Read weights of the subsequent energy terms.
106 call card_concat(weightcard)
107 call reada(weightcard,'WSC',wsc,1.0d0)
108 call reada(weightcard,'WLONG',wsc,wsc)
109 call reada(weightcard,'WSCP',wscp,1.0d0)
110 call reada(weightcard,'WELEC',welec,1.0D0)
111 call reada(weightcard,'WVDWPP',wvdwpp,welec)
112 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
113 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
114 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
115 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
116 call reada(weightcard,'WTURN3',wturn3,1.0D0)
117 call reada(weightcard,'WTURN4',wturn4,1.0D0)
118 call reada(weightcard,'WTURN6',wturn6,1.0D0)
119 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
120 call reada(weightcard,'WBOND',wbond,1.0D0)
121 call reada(weightcard,'WTOR',wtor,1.0D0)
122 call reada(weightcard,'WTORD',wtor_d,1.0D0)
123 call reada(weightcard,'WANG',wang,1.0D0)
124 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
125 call reada(weightcard,'SCAL14',scal14,0.4D0)
126 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
127 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
128 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
129 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
130 call reada(weightcard,"D0CM",d0cm,3.78d0)
131 call reada(weightcard,"AKCM",akcm,15.1d0)
132 call reada(weightcard,"AKTH",akth,11.0d0)
133 call reada(weightcard,"AKCT",akct,12.0d0)
134 call reada(weightcard,"V1SS",v1ss,-1.08d0)
135 call reada(weightcard,"V2SS",v2ss,7.61d0)
136 call reada(weightcard,"V3SS",v3ss,13.7d0)
137 call reada(weightcard,"EBR",ebr,-5.50D0)
138 if (index(weightcard,'SOFT').gt.0) ipot=6
139 C 12/1/95 Added weight for the multi-body term WCORR
140 call reada(weightcard,'WCORRH',wcorr,1.0D0)
143 dyn_ssbond_ij(i,j)=1.0d300
146 call reada(weightcard,"HT",Ht,0.0D0)
148 ss_depth=ebr/wsc-0.25*eps(1,1)
149 Ht=Ht/wsc-0.25*eps(1,1)
150 akcm=akcm*wstrain/wsc
151 akth=akth*wstrain/wsc
152 akct=akct*wstrain/wsc
153 v1ss=v1ss*wstrain/wsc
154 v2ss=v2ss*wstrain/wsc
155 v3ss=v3ss*wstrain/wsc
157 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
159 write (iout,'(/a)') "Disulfide bridge parameters:"
160 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
161 write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
162 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
163 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
164 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
166 write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
167 if (wcorr4.gt.0.0d0) wcorr=wcorr4
186 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
187 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
188 & wturn4,wturn6,wsccor
189 10 format (/'Energy-term weights (unscaled):'//
190 & 'WSCC= ',f10.6,' (SC-SC)'/
191 & 'WSCP= ',f10.6,' (SC-p)'/
192 & 'WELEC= ',f10.6,' (p-p electr)'/
193 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
194 & 'WBOND= ',f10.6,' (stretching)'/
195 & 'WANG= ',f10.6,' (bending)'/
196 & 'WSCLOC= ',f10.6,' (SC local)'/
197 & 'WTOR= ',f10.6,' (torsional)'/
198 & 'WTORD= ',f10.6,' (double torsional)'/
199 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
200 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
201 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
202 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
203 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
204 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
205 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
206 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
207 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
208 if (wcorr4.gt.0.0d0) then
209 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
210 & 'between contact pairs of peptide groups'
211 write (iout,'(2(a,f5.3/))')
212 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
213 & 'Range of quenching the correlation terms:',2*delt_corr
214 else if (wcorr.gt.0.0d0) then
215 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
216 & 'between contact pairs of peptide groups'
218 write (iout,'(a,f8.3)')
219 & 'Scaling factor of 1,4 SC-p interactions:',scal14
220 write (iout,'(a,f8.3)')
221 & 'General scaling factor of SC-p interactions:',scalscp
222 r0_corr=cutoff_corr-delt_corr
224 aad(i,1)=scalscp*aad(i,1)
225 aad(i,2)=scalscp*aad(i,2)
226 bad(i,1)=scalscp*bad(i,1)
227 bad(i,2)=scalscp*bad(i,2)
231 print *,'indpdb=',indpdb,' pdbref=',pdbref
233 C Read sequence if not taken from the pdb file.
234 if (iscode.gt.0) then
235 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
237 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
239 C Convert sequence to numeric code
241 itype(i)=rescode(i,sequence(i),iscode)
244 print '(20i4)',(itype(i),i=1,nres)
248 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
250 if (itype(i).eq.21) then
254 else if (itype(i+1).ne.20) then
256 else if (itype(i).ne.20) then
263 write (iout,*) "ITEL"
265 write (iout,*) i,itype(i),itel(i)
268 print *,'Call Read_Bridge.'
270 if (with_dihed_constr) then
272 read (inp,*) ndih_constr
273 if (ndih_constr.gt.0) then
275 write (iout,*) 'FTORS',ftors
276 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
278 & 'There are',ndih_constr,' constraints on phi angles.'
280 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
283 phi0(i)=deg2rad*phi0(i)
284 drange(i)=deg2rad*drange(i)
291 print *,'NNT=',NNT,' NCT=',NCT
292 if (itype(1).eq.21) nnt=2
293 if (itype(nres).eq.21) nct=nct-1
294 if (nstart.lt.nnt) nstart=nnt
295 if (nend.gt.nct .or. nend.eq.0) nend=nct
296 write (iout,*) "nstart",nstart," nend",nend
299 c read(inp,'(a)') pdbfile
300 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
301 c open(ipdbin,file=pdbfile,status='old',err=33)
303 c 33 write (iout,'(a)') 'Error opening PDB file.'
306 c print *,'Begin reading pdb data'
308 c print *,'Finished reading pdb data'
309 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
311 c itype_pdb(i)=itype(i)
314 c write (iout,'(a,i3)') 'nsup=',nsup
316 c if (nsup.le.(nct-nnt+1)) then
317 c do i=0,nct-nnt+1-nsup
318 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
324 c & 'Error - sequences to be superposed do not match.'
327 c do i=0,nsup-(nct-nnt+1)
328 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
330 c nstart_sup=nstart_sup+i
336 c & 'Error - sequences to be superposed do not match.'
339 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
340 c & ' nstart_seq=',nstart_seq
344 write (iout,*) "molread: REFSTR",refstr
346 if (.not.pdbref) then
347 call read_angles(inp,*38)
349 38 write (iout,'(a)') 'Error reading reference structure.'
351 call mp_stopall(Error_Msg)
353 stop 'Error reading reference structure'
365 call contact(.true.,ncont_ref,icont_ref)
367 c Read distance restraints
368 if (constr_dist.gt.0) then
369 call read_dist_constr
374 c-----------------------------------------------------------------------------
375 logical function seq_comp(itypea,itypeb,length)
377 integer length,itypea(length),itypeb(length)
380 if (itypea(i).ne.itypeb(i)) then
388 c-----------------------------------------------------------------------------
389 subroutine read_bridge
390 C Read information about disulfide bridges.
393 include 'COMMON.IOUNITS'
396 include 'COMMON.INTERACT'
397 include 'COMMON.LOCAL'
398 include 'COMMON.NAMES'
399 include 'COMMON.CHAIN'
400 include 'COMMON.FFIELD'
401 include 'COMMON.SBRIDGE'
402 include 'COMMON.HEADER'
403 include 'COMMON.CONTROL'
404 include 'COMMON.TIME1'
406 include 'COMMON.INFO'
409 C Read bridging residues.
410 read (inp,*) ns,(iss(i),i=1,ns)
412 C Check whether the specified bridging residues are cystines.
414 if (itype(iss(i)).ne.1) then
415 write (iout,'(2a,i3,a)')
416 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
417 & ' can form a disulfide bridge?!!!'
418 write (*,'(2a,i3,a)')
419 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
420 & ' can form a disulfide bridge?!!!'
422 call mp_stopall(error_msg)
428 C Read preformed bridges.
430 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
433 C Check if the residues involved in bridges are in the specified list of
437 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
438 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
439 write (iout,'(a,i3,a)') 'Disulfide pair',i,
440 & ' contains residues present in other pairs.'
441 write (*,'(a,i3,a)') 'Disulfide pair',i,
442 & ' contains residues present in other pairs.'
444 call mp_stopall(error_msg)
451 if (ihpb(i).eq.iss(j)) goto 10
453 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
456 if (jhpb(i).eq.iss(j)) goto 20
458 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
469 if (ns.gt.0.and.dyn_ss) then
473 forcon(i-nss)=forcon(i)
480 dyn_ss_mask(iss(i))=.true.
481 c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
484 print *, "Leaving brigde read"
487 c----------------------------------------------------------------------------
488 subroutine read_angles(kanal,*)
493 include 'COMMON.CHAIN'
494 include 'COMMON.IOUNITS'
496 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
497 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
498 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
499 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
501 theta(i)=deg2rad*theta(i)
502 phi(i)=deg2rad*phi(i)
503 alph(i)=deg2rad*alph(i)
504 omeg(i)=deg2rad*omeg(i)
509 c----------------------------------------------------------------------------
510 subroutine reada(rekord,lancuch,wartosc,default)
512 character*(*) rekord,lancuch
513 double precision wartosc,default
516 iread=index(rekord,lancuch)
521 iread=iread+ilen(lancuch)+1
522 read (rekord(iread:),*) wartosc
525 c----------------------------------------------------------------------------
526 subroutine multreada(rekord,lancuch,tablica,dim,default)
529 double precision tablica(dim),default
530 character*(*) rekord,lancuch
536 iread=index(rekord,lancuch)
537 if (iread.eq.0) return
538 iread=iread+ilen(lancuch)+1
539 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
542 c----------------------------------------------------------------------------
543 subroutine readi(rekord,lancuch,wartosc,default)
545 character*(*) rekord,lancuch
546 integer wartosc,default
549 iread=index(rekord,lancuch)
554 iread=iread+ilen(lancuch)+1
555 read (rekord(iread:),*) wartosc
558 c----------------------------------------------------------------------------
559 subroutine multreadi(rekord,lancuch,tablica,dim,default)
562 integer tablica(dim),default
563 character*(*) rekord,lancuch
570 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
571 if (iread.eq.0) return
572 iread=iread+ilen(lancuch)+1
573 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
576 c----------------------------------------------------------------------------
577 subroutine card_concat(card)
579 include 'COMMON.IOUNITS'
581 character*80 karta,ucase
583 read (inp,'(a)') karta
586 do while (karta(80:80).eq.'&')
587 card=card(:ilen(card)+1)//karta(:79)
588 read (inp,'(a)') karta
591 card=card(:ilen(card)+1)//karta
594 c----------------------------------------------------------------------------
603 include 'COMMON.IOUNITS'
604 include 'COMMON.CONTROL'
605 integer lenpre,lenpot,ilen
607 character*16 cformat,cprint
609 integer lenint,lenout
610 call getenv('INPUT',prefix)
611 call getenv('OUTPUT',prefout)
612 call getenv('INTIN',prefintin)
613 call getenv('COORD',cformat)
614 call getenv('PRINTCOOR',cprint)
615 call getenv('SCRATCHDIR',scratchdir)
618 if (index(ucase(cformat),'CX').gt.0) then
625 lenint=ilen(prefintin)
626 C Get the names and open the input files
627 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
629 write (liczba,'(bz,i3.3)') me
630 outname=prefout(:lenout)//'_clust.out_'//liczba
632 outname=prefout(:lenout)//'_clust.out'
635 intinname=prefintin(:lenint)//'.bx'
636 else if (from_cx) then
637 intinname=prefintin(:lenint)//'.cx'
639 intinname=prefintin(:lenint)//'.int'
641 rmsname=prefintin(:lenint)//'.rms'
642 open (jplot,file=prefout(:ilen(prefout))//'.tex',
644 open (jrms,file=rmsname,status='unknown')
645 open(iout,file=outname,status='unknown')
646 C Get parameter filenames and open the parameter files.
647 call getenv('BONDPAR',bondname)
648 open (ibond,file=bondname,status='old')
649 call getenv('THETPAR',thetname)
650 open (ithep,file=thetname,status='old')
651 call getenv('ROTPAR',rotname)
652 open (irotam,file=rotname,status='old')
653 call getenv('TORPAR',torname)
654 open (itorp,file=torname,status='old')
655 call getenv('TORDPAR',tordname)
656 open (itordp,file=tordname,status='old')
657 call getenv('FOURIER',fouriername)
658 open (ifourier,file=fouriername,status='old')
659 call getenv('ELEPAR',elename)
660 open (ielep,file=elename,status='old')
661 call getenv('SIDEPAR',sidename)
662 open (isidep,file=sidename,status='old')
663 call getenv('SIDEP',sidepname)
664 open (isidep1,file=sidepname,status="old")
665 call getenv('SCCORPAR',sccorname)
666 open (isccor,file=sccorname,status="old")
669 C 8/9/01 In the newest version SCp interaction constants are read from a file
670 C Use -DOLDSCP to use hard-coded constants instead.
672 call getenv('SCPPAR',scpname)
673 open (iscpp,file=scpname,status='old')
677 c-------------------------------------------------------------------------------
678 subroutine read_dist_constr
679 implicit real*8 (a-h,o-z)
681 include 'COMMON.CONTROL'
682 include 'COMMON.CHAIN'
683 include 'COMMON.IOUNITS'
684 include 'COMMON.SBRIDGE'
685 integer ifrag_(2,100),ipair_(2,100)
686 double precision wfrag_(100),wpair_(100)
687 character*500 controlcard
688 c write (iout,*) "Calling read_dist_constr"
689 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
690 call card_concat(controlcard)
692 c write (iout,'(a)') controlcard
693 call readi(controlcard,"NFRAG",nfrag_,0)
694 call readi(controlcard,"NPAIR",npair_,0)
695 call readi(controlcard,"NDIST",ndist_,0)
696 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
697 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
698 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
699 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
700 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
701 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
702 write (iout,*) "IFRAG"
704 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
706 write (iout,*) "IPAIR"
708 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
711 if (.not.refstr .and. nfrag_.gt.0) then
713 & "ERROR: no reference structure to compute distance restraints"
715 & "Restraints must be specified explicitly (NDIST=number)"
718 if (nfrag_.lt.2 .and. npair_.gt.0) then
719 write (iout,*) "ERROR: Less than 2 fragments specified",
720 & " but distance restraints between pairs requested"
725 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
726 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
727 & ifrag_(2,i)=nstart_sup+nsup-1
728 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
730 if (wfrag_(i).gt.0.0d0) then
731 do j=ifrag_(1,i),ifrag_(2,i)-1
733 write (iout,*) "j",j," k",k
735 if (constr_dist.eq.1) then
740 forcon(nhpb)=wfrag_(i)
741 else if (constr_dist.eq.2) then
742 if (ddjk.le.dist_cut) then
747 forcon(nhpb)=wfrag_(i)
754 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
756 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
757 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
763 if (wpair_(i).gt.0.0d0) then
771 do j=ifrag_(1,ii),ifrag_(2,ii)
772 do k=ifrag_(1,jj),ifrag_(2,jj)
776 forcon(nhpb)=wpair_(i)
778 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
779 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
785 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
786 & ibecarb(i),forcon(nhpb+1)
787 if (forcon(nhpb+1).gt.0.0d0) then
789 if (ibecarb(i).gt.0) then
793 if (dhpb(nhpb).eq.0.0d0)
794 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
798 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
799 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)