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 character*320 controlcard,ucase
22 integer i,i1,i2,it1,it2
24 read (INP,'(a80)') titel
25 call card_concat(controlcard)
27 call readi(controlcard,'NRES',nres,0)
28 call readi(controlcard,'RESCALE',rescale_mode,2)
29 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
30 write (iout,*) "DISTCHAINMAX",distchainmax
31 call readi(controlcard,'PDBOUT',outpdb,0)
32 call readi(controlcard,'MOL2OUT',outmol2,0)
33 refstr=(index(controlcard,'REFSTR').gt.0)
34 write (iout,*) "REFSTR",refstr
35 pdbref=(index(controlcard,'PDBREF').gt.0)
36 iscode=index(controlcard,'ONE_LETTER')
37 tree=(index(controlcard,'MAKE_TREE').gt.0)
38 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
39 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
40 write (iout,*) "with_dihed_constr ",with_dihed_constr,
41 & " CONSTR_DIST",constr_dist
42 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
43 write (iout,*) "with_theta_constr ",with_theta_constr
45 min_var=(index(controlcard,'MINVAR').gt.0)
46 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
47 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
48 call readi(controlcard,'NCUT',ncut,1)
49 call readi(controlcard,'SYM',symetr,1)
50 write (iout,*) 'sym', symetr
51 call readi(controlcard,'NSTART',nstart,0)
52 call readi(controlcard,'NEND',nend,0)
53 call reada(controlcard,'ECUT',ecut,10.0d0)
54 call reada(controlcard,'PROB',prob_limit,0.99d0)
55 write (iout,*) "Probability limit",prob_limit
56 lgrp=(index(controlcard,'LGRP').gt.0)
57 caonly=(index(controlcard,'CA_ONLY').gt.0)
58 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
59 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
60 call readi(controlcard,'IOPT',iopt,2)
61 lside = index(controlcard,"SIDE").gt.0
62 efree = index(controlcard,"EFREE").gt.0
63 call readi(controlcard,'NTEMP',nT,1)
64 write (iout,*) "nT",nT
65 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
66 write (iout,*) "nT",nT
67 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
69 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
71 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
72 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
73 lprint_int=index(controlcard,"PRINT_INT") .gt.0
77 c--------------------------------------------------------------------------
80 C Read molecular data.
84 include 'COMMON.IOUNITS'
87 include 'COMMON.INTERACT'
88 include 'COMMON.LOCAL'
89 include 'COMMON.NAMES'
90 include 'COMMON.CHAIN'
91 include 'COMMON.FFIELD'
92 include 'COMMON.SBRIDGE'
93 include 'COMMON.HEADER'
94 include 'COMMON.CONTROL'
95 include 'COMMON.CONTACTS'
96 include 'COMMON.TIME1'
97 include 'COMMON.TORCNSTR'
101 character*4 sequence(maxres)
102 character*800 weightcard
104 double precision x(maxvar)
105 integer itype_pdb(maxres)
107 integer i,j,kkk,i1,i2,it1,it2
111 C Read weights of the subsequent energy terms.
112 call card_concat(weightcard)
113 call reada(weightcard,'WSC',wsc,1.0d0)
114 call reada(weightcard,'WLONG',wsc,wsc)
115 call reada(weightcard,'WSCP',wscp,1.0d0)
116 call reada(weightcard,'WELEC',welec,1.0D0)
117 call reada(weightcard,'WVDWPP',wvdwpp,welec)
118 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
119 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
120 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
121 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
122 call reada(weightcard,'WTURN3',wturn3,1.0D0)
123 call reada(weightcard,'WTURN4',wturn4,1.0D0)
124 call reada(weightcard,'WTURN6',wturn6,1.0D0)
125 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
126 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
127 call reada(weightcard,'WBOND',wbond,1.0D0)
128 call reada(weightcard,'WTOR',wtor,1.0D0)
129 call reada(weightcard,'WTORD',wtor_d,1.0D0)
130 call reada(weightcard,'WANG',wang,1.0D0)
131 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
132 call reada(weightcard,'SCAL14',scal14,0.4D0)
133 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
134 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
135 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
136 if (index(weightcard,'SOFT').gt.0) ipot=6
137 call reada(weightcard,"D0CM",d0cm,3.78d0)
138 call reada(weightcard,"AKCM",akcm,15.1d0)
139 call reada(weightcard,"AKTH",akth,11.0d0)
140 call reada(weightcard,"AKCT",akct,12.0d0)
141 call reada(weightcard,"V1SS",v1ss,-1.08d0)
142 call reada(weightcard,"V2SS",v2ss,7.61d0)
143 call reada(weightcard,"V3SS",v3ss,13.7d0)
144 call reada(weightcard,"EBR",ebr,-5.50D0)
145 call reada(weightcard,"ATRISS",atriss,0.301D0)
146 call reada(weightcard,"BTRISS",btriss,0.021D0)
147 call reada(weightcard,"CTRISS",ctriss,1.001D0)
148 call reada(weightcard,"DTRISS",dtriss,1.001D0)
149 write (iout,*) "ATRISS=", atriss
150 write (iout,*) "BTRISS=", btriss
151 write (iout,*) "CTRISS=", ctriss
152 write (iout,*) "DTRISS=", dtriss
153 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
155 dyn_ss_mask(i)=.false.
159 dyn_ssbond_ij(i,j)=1.0d300
162 call reada(weightcard,"HT",Ht,0.0D0)
164 ss_depth=ebr/wsc-0.25*eps(1,1)
165 Ht=Ht/wsc-0.25*eps(1,1)
166 akcm=akcm*wstrain/wsc
167 akth=akth*wstrain/wsc
168 akct=akct*wstrain/wsc
169 v1ss=v1ss*wstrain/wsc
170 v2ss=v2ss*wstrain/wsc
171 v3ss=v3ss*wstrain/wsc
173 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
175 write (iout,'(/a)') "Disulfide bridge parameters:"
176 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
177 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
178 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
179 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
182 C 12/1/95 Added weight for the multi-body term WCORR
183 call reada(weightcard,'WCORRH',wcorr,1.0D0)
184 if (wcorr4.gt.0.0d0) wcorr=wcorr4
204 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
205 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
206 & wturn4,wturn6,wsccor
207 10 format (/'Energy-term weights (unscaled):'//
208 & 'WSCC= ',f10.6,' (SC-SC)'/
209 & 'WSCP= ',f10.6,' (SC-p)'/
210 & 'WELEC= ',f10.6,' (p-p electr)'/
211 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
212 & 'WBOND= ',f10.6,' (stretching)'/
213 & 'WANG= ',f10.6,' (bending)'/
214 & 'WSCLOC= ',f10.6,' (SC local)'/
215 & 'WTOR= ',f10.6,' (torsional)'/
216 & 'WTORD= ',f10.6,' (double torsional)'/
217 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
218 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
219 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
220 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
221 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
222 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
223 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
224 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
225 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
227 if (wcorr4.gt.0.0d0) then
228 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
229 & 'between contact pairs of peptide groups'
230 write (iout,'(2(a,f5.3/))')
231 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
232 & 'Range of quenching the correlation terms:',2*delt_corr
233 else if (wcorr.gt.0.0d0) then
234 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
235 & 'between contact pairs of peptide groups'
237 write (iout,'(a,f8.3)')
238 & 'Scaling factor of 1,4 SC-p interactions:',scal14
239 write (iout,'(a,f8.3)')
240 & 'General scaling factor of SC-p interactions:',scalscp
241 r0_corr=cutoff_corr-delt_corr
243 aad(i,1)=scalscp*aad(i,1)
244 aad(i,2)=scalscp*aad(i,2)
245 bad(i,1)=scalscp*bad(i,1)
246 bad(i,2)=scalscp*bad(i,2)
250 print *,'indpdb=',indpdb,' pdbref=',pdbref
252 C Read sequence if not taken from the pdb file.
253 if (iscode.gt.0) then
254 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
256 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
258 C Convert sequence to numeric code
260 itype(i)=rescode(i,sequence(i),iscode)
263 print '(20i4)',(itype(i),i=1,nres)
267 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
269 if (itype(i).eq.ntyp1) then
273 else if (iabs(itype(i+1)).ne.20) then
275 else if (iabs(itype(i)).ne.20) then
282 write (iout,*) "ITEL"
284 write (iout,*) i,itype(i),itel(i)
287 print *,'Call Read_Bridge.'
289 C this fragment reads diheadral constrains
290 if (with_dihed_constr) then
292 read (inp,*) ndih_constr
293 if (ndih_constr.gt.0) then
295 C write (iout,*) 'FTORS',ftors
296 C ftors is the force constant for torsional quartic constrains
297 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
300 & 'There are',ndih_constr,' constraints on phi angles.'
302 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
306 phi0(i)=deg2rad*phi0(i)
307 drange(i)=deg2rad*drange(i)
309 endif ! endif ndif_constr.gt.0
310 endif ! with_dihed_constr
311 if (with_theta_constr) then
312 C with_theta_constr is keyword allowing for occurance of theta constrains
313 read (inp,*) ntheta_constr
314 C ntheta_constr is the number of theta constrains
315 if (ntheta_constr.gt.0) then
317 read (inp,*) (itheta_constr(i),theta_constr0(i),
318 & theta_drange(i),for_thet_constr(i),
320 C the above code reads from 1 to ntheta_constr
321 C itheta_constr(i) residue i for which is theta_constr
322 C theta_constr0 the global minimum value
323 C theta_drange is range for which there is no energy penalty
324 C for_thet_constr is the force constant for quartic energy penalty
326 C if(me.eq.king.or..not.out1file)then
328 & 'There are',ntheta_constr,' constraints on phi angles.'
330 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
336 theta_constr0(i)=deg2rad*theta_constr0(i)
337 theta_drange(i)=deg2rad*theta_drange(i)
339 C if(me.eq.king.or..not.out1file)
340 C & write (iout,*) 'FTORS',ftors
341 C do i=1,ntheta_constr
342 C ii = itheta_constr(i)
343 C thetabound(1,ii) = phi0(i)-drange(i)
344 C thetabound(2,ii) = phi0(i)+drange(i)
346 endif ! ntheta_constr.gt.0
347 endif! with_theta_constr
351 print *,'NNT=',NNT,' NCT=',NCT
352 if (itype(1).eq.ntyp1) nnt=2
353 if (itype(nres).eq.ntyp1) nct=nct-1
354 if (nstart.lt.nnt) nstart=nnt
355 if (nend.gt.nct .or. nend.eq.0) nend=nct
356 write (iout,*) "nstart",nstart," nend",nend
359 c read(inp,'(a)') pdbfile
360 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
361 c open(ipdbin,file=pdbfile,status='old',err=33)
363 c 33 write (iout,'(a)') 'Error opening PDB file.'
366 c print *,'Begin reading pdb data'
368 c print *,'Finished reading pdb data'
369 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
371 c itype_pdb(i)=itype(i)
374 c write (iout,'(a,i3)') 'nsup=',nsup
376 c if (nsup.le.(nct-nnt+1)) then
377 c do i=0,nct-nnt+1-nsup
378 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
384 c & 'Error - sequences to be superposed do not match.'
387 c do i=0,nsup-(nct-nnt+1)
388 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
390 c nstart_sup=nstart_sup+i
396 c & 'Error - sequences to be superposed do not match.'
399 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
400 c & ' nstart_seq=',nstart_seq
404 write (iout,*) "molread: REFSTR",refstr
406 if (.not.pdbref) then
407 call read_angles(inp,*38)
409 38 write (iout,'(a)') 'Error reading reference structure.'
411 call mp_stopall(Error_Msg)
413 stop 'Error reading reference structure'
426 call contact(.true.,ncont_ref,icont_ref)
429 C write (iout,'(/a,i3,a)')
430 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
431 write (iout,'(20i4)') (iss(i),i=1,ns)
433 write(iout,*)"Running with dynamic disulfide-bond formation"
435 write (iout,'(/a/)') 'Pre-formed links are:'
441 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
442 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
448 if (ns.gt.0.and.dyn_ss) then
452 forcon(i-nss)=forcon(i)
459 dyn_ss_mask(iss(i))=.true.
462 c Read distance restraints
463 if (constr_dist.gt.0) then
464 call read_dist_constr
469 c-----------------------------------------------------------------------------
470 logical function seq_comp(itypea,itypeb,length)
472 integer length,itypea(length),itypeb(length)
475 if (itypea(i).ne.itypeb(i)) then
483 c-----------------------------------------------------------------------------
484 subroutine read_bridge
485 C Read information about disulfide bridges.
488 include 'COMMON.IOUNITS'
491 include 'COMMON.INTERACT'
492 include 'COMMON.LOCAL'
493 include 'COMMON.NAMES'
494 include 'COMMON.CHAIN'
495 include 'COMMON.FFIELD'
496 include 'COMMON.SBRIDGE'
497 include 'COMMON.HEADER'
498 include 'COMMON.CONTROL'
499 include 'COMMON.TIME1'
501 include 'COMMON.INFO'
504 C Read bridging residues.
505 read (inp,*) ns,(iss(i),i=1,ns)
507 C Check whether the specified bridging residues are cystines.
509 if (itype(iss(i)).ne.1) then
510 write (iout,'(2a,i3,a)')
511 & 'Do you REALLY think that the residue ',
512 & restyp(itype(iss(i))),i,
513 & ' can form a disulfide bridge?!!!'
514 write (*,'(2a,i3,a)')
515 & 'Do you REALLY think that the residue ',
516 & restyp(itype(iss(i))),i,
517 & ' can form a disulfide bridge?!!!'
519 call mp_stopall(error_msg)
525 C Read preformed bridges.
527 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
530 C Check if the residues involved in bridges are in the specified list of
534 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
535 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
536 write (iout,'(a,i3,a)') 'Disulfide pair',i,
537 & ' contains residues present in other pairs.'
538 write (*,'(a,i3,a)') 'Disulfide pair',i,
539 & ' contains residues present in other pairs.'
541 call mp_stopall(error_msg)
548 if (ihpb(i).eq.iss(j)) goto 10
550 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
553 if (jhpb(i).eq.iss(j)) goto 20
555 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
568 c----------------------------------------------------------------------------
569 subroutine read_angles(kanal,*)
574 include 'COMMON.CHAIN'
575 include 'COMMON.IOUNITS'
577 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
578 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
579 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
580 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
582 theta(i)=deg2rad*theta(i)
583 phi(i)=deg2rad*phi(i)
584 alph(i)=deg2rad*alph(i)
585 omeg(i)=deg2rad*omeg(i)
590 c----------------------------------------------------------------------------
591 subroutine reada(rekord,lancuch,wartosc,default)
593 character*(*) rekord,lancuch
594 double precision wartosc,default
597 iread=index(rekord,lancuch)
602 iread=iread+ilen(lancuch)+1
603 read (rekord(iread:),*) wartosc
606 c----------------------------------------------------------------------------
607 subroutine multreada(rekord,lancuch,tablica,dim,default)
610 double precision tablica(dim),default
611 character*(*) rekord,lancuch
617 iread=index(rekord,lancuch)
618 if (iread.eq.0) return
619 iread=iread+ilen(lancuch)+1
620 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
623 c----------------------------------------------------------------------------
624 subroutine readi(rekord,lancuch,wartosc,default)
626 character*(*) rekord,lancuch
627 integer wartosc,default
630 iread=index(rekord,lancuch)
635 iread=iread+ilen(lancuch)+1
636 read (rekord(iread:),*) wartosc
639 C----------------------------------------------------------------------
640 subroutine multreadi(rekord,lancuch,tablica,dim,default)
643 integer tablica(dim),default
644 character*(*) rekord,lancuch
651 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
652 if (iread.eq.0) return
653 iread=iread+ilen(lancuch)+1
654 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
658 c----------------------------------------------------------------------------
659 subroutine card_concat(card)
661 include 'COMMON.IOUNITS'
663 character*80 karta,ucase
665 read (inp,'(a)') karta
668 do while (karta(80:80).eq.'&')
669 card=card(:ilen(card)+1)//karta(:79)
670 read (inp,'(a)') karta
673 card=card(:ilen(card)+1)//karta
676 c----------------------------------------------------------------------------
685 include 'COMMON.IOUNITS'
686 include 'COMMON.CONTROL'
687 integer lenpre,lenpot,ilen
689 character*16 cformat,cprint
691 integer lenint,lenout
692 call getenv('INPUT',prefix)
693 call getenv('OUTPUT',prefout)
694 call getenv('INTIN',prefintin)
695 call getenv('COORD',cformat)
696 call getenv('PRINTCOOR',cprint)
697 call getenv('SCRATCHDIR',scratchdir)
700 if (index(ucase(cformat),'CX').gt.0) then
707 lenint=ilen(prefintin)
708 C Get the names and open the input files
709 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
711 write (liczba,'(bz,i3.3)') me
712 outname=prefout(:lenout)//'_clust.out_'//liczba
714 outname=prefout(:lenout)//'_clust.out'
717 intinname=prefintin(:lenint)//'.bx'
718 else if (from_cx) then
719 intinname=prefintin(:lenint)//'.cx'
721 intinname=prefintin(:lenint)//'.int'
723 rmsname=prefintin(:lenint)//'.rms'
724 open (jplot,file=prefout(:ilen(prefout))//'.tex',
726 open (jrms,file=rmsname,status='unknown')
727 open(iout,file=outname,status='unknown')
728 C Get parameter filenames and open the parameter files.
729 call getenv('BONDPAR',bondname)
730 open (ibond,file=bondname,status='old')
731 call getenv('THETPAR',thetname)
732 open (ithep,file=thetname,status='old')
733 call getenv('ROTPAR',rotname)
734 open (irotam,file=rotname,status='old')
735 call getenv('TORPAR',torname)
736 open (itorp,file=torname,status='old')
737 call getenv('TORDPAR',tordname)
738 open (itordp,file=tordname,status='old')
739 call getenv('FOURIER',fouriername)
740 open (ifourier,file=fouriername,status='old')
741 call getenv('ELEPAR',elename)
742 open (ielep,file=elename,status='old')
743 call getenv('SIDEPAR',sidename)
744 open (isidep,file=sidename,status='old')
745 call getenv('SIDEP',sidepname)
746 open (isidep1,file=sidepname,status="old")
747 call getenv('SCCORPAR',sccorname)
748 open (isccor,file=sccorname,status="old")
751 C 8/9/01 In the newest version SCp interaction constants are read from a file
752 C Use -DOLDSCP to use hard-coded constants instead.
754 call getenv('SCPPAR',scpname)
755 open (iscpp,file=scpname,status='old')
759 subroutine read_dist_constr
760 implicit real*8 (a-h,o-z)
765 include 'COMMON.CONTROL'
766 include 'COMMON.CHAIN'
767 include 'COMMON.IOUNITS'
768 include 'COMMON.SBRIDGE'
769 integer ifrag_(2,100),ipair_(2,100)
770 double precision wfrag_(100),wpair_(100)
771 character*500 controlcard
772 logical lprn /.true./
773 write (iout,*) "Calling read_dist_constr"
774 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
776 write(iout,*) "TU sie wywalam?"
777 call card_concat(controlcard)
778 write (iout,*) controlcard
780 call readi(controlcard,"NFRAG",nfrag_,0)
781 call readi(controlcard,"NPAIR",npair_,0)
782 call readi(controlcard,"NDIST",ndist_,0)
783 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
784 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
785 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
786 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
787 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
788 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
789 write (iout,*) "IFRAG"
791 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
793 write (iout,*) "IPAIR"
795 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
799 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
800 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
801 & ifrag_(2,i)=nstart_sup+nsup-1
802 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
804 if (wfrag_(i).gt.0.0d0) then
805 do j=ifrag_(1,i),ifrag_(2,i)-1
807 write (iout,*) "j",j," k",k
809 if (constr_dist.eq.1) then
814 forcon(nhpb)=wfrag_(i)
815 else if (constr_dist.eq.2) then
816 if (ddjk.le.dist_cut) then
821 forcon(nhpb)=wfrag_(i)
828 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
831 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
832 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
838 if (wpair_(i).gt.0.0d0) then
846 do j=ifrag_(1,ii),ifrag_(2,ii)
847 do k=ifrag_(1,jj),ifrag_(2,jj)
851 forcon(nhpb)=wpair_(i)
853 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
854 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
860 if (constr_dist.eq.11) then
861 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
862 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
863 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
864 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
865 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
867 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
869 if (forcon(nhpb+1).gt.0.0d0) then
871 if (ibecarb(i).gt.0) then
875 if (dhpb(nhpb).eq.0.0d0)
876 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
877 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
878 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
879 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)