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 include 'COMMON.SHIELD'
20 character*320 controlcard,ucase
24 integer i,i1,i2,it1,it2
26 read (INP,'(a80)') titel
27 call card_concat(controlcard)
29 call readi(controlcard,'NRES',nres,0)
30 call readi(controlcard,'RESCALE',rescale_mode,2)
31 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
32 write (iout,*) "DISTCHAINMAX",distchainmax
33 C Reading the dimensions of box in x,y,z coordinates
34 call reada(controlcard,'BOXX',boxxsize,100.0d0)
35 call reada(controlcard,'BOXY',boxysize,100.0d0)
36 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
37 c Cutoff range for interactions
38 call reada(controlcard,"R_CUT",r_cut,15.0d0)
39 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
41 call readi(controlcard,'SHIELD',shield_mode,0)
42 write (iout,*) "SHIELD MODE",shield_mode
43 if (shield_mode.gt.0) then
45 C VSolvSphere the volume of solving sphere
47 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
48 C there will be no distinction between proline peptide group and normal peptide
49 C group in case of shielding parameters
50 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
51 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
52 write (iout,*) VSolvSphere,VSolvSphere_div
53 C long axis of side chain
55 long_r_sidechain(i)=vbldsc0(1,i)
56 short_r_sidechain(i)=sigma0(i)
60 call readi(controlcard,'PDBOUT',outpdb,0)
61 call readi(controlcard,'MOL2OUT',outmol2,0)
62 refstr=(index(controlcard,'REFSTR').gt.0)
63 write (iout,*) "REFSTR",refstr
64 pdbref=(index(controlcard,'PDBREF').gt.0)
65 iscode=index(controlcard,'ONE_LETTER')
66 tree=(index(controlcard,'MAKE_TREE').gt.0)
67 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
68 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
69 write (iout,*) "with_dihed_constr ",with_dihed_constr,
70 & " CONSTR_DIST",constr_dist
71 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
72 write (iout,*) "with_theta_constr ",with_theta_constr
74 min_var=(index(controlcard,'MINVAR').gt.0)
75 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
76 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
77 call readi(controlcard,'NCUT',ncut,1)
78 call readi(controlcard,'SYM',symetr,1)
79 write (iout,*) 'sym', symetr
80 call readi(controlcard,'NSTART',nstart,0)
81 call readi(controlcard,'NEND',nend,0)
82 call reada(controlcard,'ECUT',ecut,10.0d0)
83 call reada(controlcard,'PROB',prob_limit,0.99d0)
84 write (iout,*) "Probability limit",prob_limit
85 lgrp=(index(controlcard,'LGRP').gt.0)
86 caonly=(index(controlcard,'CA_ONLY').gt.0)
87 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
88 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
89 call readi(controlcard,'IOPT',iopt,2)
90 lside = index(controlcard,"SIDE").gt.0
91 efree = index(controlcard,"EFREE").gt.0
92 call readi(controlcard,'NTEMP',nT,1)
93 write (iout,*) "nT",nT
94 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
95 write (iout,*) "nT",nT
96 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
98 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
100 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
101 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
102 lprint_int=index(controlcard,"PRINT_INT") .gt.0
106 c--------------------------------------------------------------------------
109 C Read molecular data.
113 include 'COMMON.IOUNITS'
116 include 'COMMON.INTERACT'
117 include 'COMMON.LOCAL'
118 include 'COMMON.NAMES'
119 include 'COMMON.CHAIN'
120 include 'COMMON.FFIELD'
121 include 'COMMON.SBRIDGE'
122 include 'COMMON.HEADER'
123 include 'COMMON.CONTROL'
124 include 'COMMON.CONTACTS'
125 include 'COMMON.TIME1'
126 include 'COMMON.TORCNSTR'
127 include 'COMMON.SHIELD'
129 include 'COMMON.INFO'
131 character*4 sequence(maxres)
132 character*800 weightcard
134 double precision x(maxvar)
135 integer itype_pdb(maxres)
137 integer i,j,kkk,i1,i2,it1,it2
141 C Read weights of the subsequent energy terms.
142 call card_concat(weightcard)
143 write(iout,*) weightcard
144 C call reada(weightcard,'WSC',wsc,1.0d0)
146 call reada(weightcard,'WLONG',wsc,wsc)
147 call reada(weightcard,'WSCP',wscp,1.0d0)
148 call reada(weightcard,'WELEC',welec,1.0D0)
149 call reada(weightcard,'WVDWPP',wvdwpp,welec)
150 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
151 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
152 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
153 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
154 call reada(weightcard,'WTURN3',wturn3,1.0D0)
155 call reada(weightcard,'WTURN4',wturn4,1.0D0)
156 call reada(weightcard,'WTURN6',wturn6,1.0D0)
157 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
158 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
159 call reada(weightcard,'WBOND',wbond,1.0D0)
160 call reada(weightcard,'WTOR',wtor,1.0D0)
161 call reada(weightcard,'WTORD',wtor_d,1.0D0)
162 call reada(weightcard,'WANG',wang,1.0D0)
163 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
164 call reada(weightcard,'SCAL14',scal14,0.4D0)
165 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
166 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
167 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
168 if (index(weightcard,'SOFT').gt.0) ipot=6
169 call reada(weightcard,"D0CM",d0cm,3.78d0)
170 call reada(weightcard,"AKCM",akcm,15.1d0)
171 call reada(weightcard,"AKTH",akth,11.0d0)
172 call reada(weightcard,"AKCT",akct,12.0d0)
173 call reada(weightcard,"V1SS",v1ss,-1.08d0)
174 call reada(weightcard,"V2SS",v2ss,7.61d0)
175 call reada(weightcard,"V3SS",v3ss,13.7d0)
176 call reada(weightcard,"EBR",ebr,-5.50D0)
177 call reada(weightcard,'WSHIELD',wshield,1.0d0)
178 write(iout,*) 'WSHIELD',wshield
179 call reada(weightcard,"ATRISS",atriss,0.301D0)
180 call reada(weightcard,"BTRISS",btriss,0.021D0)
181 call reada(weightcard,"CTRISS",ctriss,1.001D0)
182 call reada(weightcard,"DTRISS",dtriss,1.001D0)
183 write (iout,*) "ATRISS=", atriss
184 write (iout,*) "BTRISS=", btriss
185 write (iout,*) "CTRISS=", ctriss
186 write (iout,*) "DTRISS=", dtriss
187 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
189 dyn_ss_mask(i)=.false.
193 dyn_ssbond_ij(i,j)=1.0d300
196 call reada(weightcard,"HT",Ht,0.0D0)
198 ss_depth=ebr/wsc-0.25*eps(1,1)
199 Ht=Ht/wsc-0.25*eps(1,1)
200 akcm=akcm*wstrain/wsc
201 akth=akth*wstrain/wsc
202 akct=akct*wstrain/wsc
203 v1ss=v1ss*wstrain/wsc
204 v2ss=v2ss*wstrain/wsc
205 v3ss=v3ss*wstrain/wsc
207 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
209 write (iout,'(/a)') "Disulfide bridge parameters:"
210 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
211 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
212 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
213 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
216 C 12/1/95 Added weight for the multi-body term WCORR
217 call reada(weightcard,'WCORRH',wcorr,1.0D0)
218 if (wcorr4.gt.0.0d0) wcorr=wcorr4
238 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
239 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
240 & wturn4,wturn6,wsccor
241 10 format (/'Energy-term weights (unscaled):'//
242 & 'WSCC= ',f10.6,' (SC-SC)'/
243 & 'WSCP= ',f10.6,' (SC-p)'/
244 & 'WELEC= ',f10.6,' (p-p electr)'/
245 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
246 & 'WBOND= ',f10.6,' (stretching)'/
247 & 'WANG= ',f10.6,' (bending)'/
248 & 'WSCLOC= ',f10.6,' (SC local)'/
249 & 'WTOR= ',f10.6,' (torsional)'/
250 & 'WTORD= ',f10.6,' (double torsional)'/
251 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
252 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
253 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
254 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
255 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
256 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
257 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
258 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
259 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
261 if (wcorr4.gt.0.0d0) then
262 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
263 & 'between contact pairs of peptide groups'
264 write (iout,'(2(a,f5.3/))')
265 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
266 & 'Range of quenching the correlation terms:',2*delt_corr
267 else if (wcorr.gt.0.0d0) then
268 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
269 & 'between contact pairs of peptide groups'
271 write (iout,'(a,f8.3)')
272 & 'Scaling factor of 1,4 SC-p interactions:',scal14
273 write (iout,'(a,f8.3)')
274 & 'General scaling factor of SC-p interactions:',scalscp
275 r0_corr=cutoff_corr-delt_corr
277 aad(i,1)=scalscp*aad(i,1)
278 aad(i,2)=scalscp*aad(i,2)
279 bad(i,1)=scalscp*bad(i,1)
280 bad(i,2)=scalscp*bad(i,2)
284 print *,'indpdb=',indpdb,' pdbref=',pdbref
286 C Read sequence if not taken from the pdb file.
287 if (iscode.gt.0) then
288 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
290 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
292 C Convert sequence to numeric code
294 itype(i)=rescode(i,sequence(i),iscode)
297 print '(20i4)',(itype(i),i=1,nres)
301 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
303 if (itype(i).eq.ntyp1) then
307 else if (iabs(itype(i+1)).ne.20) then
309 else if (iabs(itype(i)).ne.20) then
316 write (iout,*) "ITEL"
318 write (iout,*) i,itype(i),itel(i)
321 print *,'Call Read_Bridge.'
323 C this fragment reads diheadral constrains
324 if (with_dihed_constr) then
326 read (inp,*) ndih_constr
327 if (ndih_constr.gt.0) then
329 C write (iout,*) 'FTORS',ftors
330 C ftors is the force constant for torsional quartic constrains
331 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
334 & 'There are',ndih_constr,' constraints on phi angles.'
336 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
340 phi0(i)=deg2rad*phi0(i)
341 drange(i)=deg2rad*drange(i)
343 endif ! endif ndif_constr.gt.0
344 endif ! with_dihed_constr
345 if (with_theta_constr) then
346 C with_theta_constr is keyword allowing for occurance of theta constrains
347 read (inp,*) ntheta_constr
348 C ntheta_constr is the number of theta constrains
349 if (ntheta_constr.gt.0) then
351 read (inp,*) (itheta_constr(i),theta_constr0(i),
352 & theta_drange(i),for_thet_constr(i),
354 C the above code reads from 1 to ntheta_constr
355 C itheta_constr(i) residue i for which is theta_constr
356 C theta_constr0 the global minimum value
357 C theta_drange is range for which there is no energy penalty
358 C for_thet_constr is the force constant for quartic energy penalty
360 C if(me.eq.king.or..not.out1file)then
362 & 'There are',ntheta_constr,' constraints on phi angles.'
364 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
370 theta_constr0(i)=deg2rad*theta_constr0(i)
371 theta_drange(i)=deg2rad*theta_drange(i)
373 C if(me.eq.king.or..not.out1file)
374 C & write (iout,*) 'FTORS',ftors
375 C do i=1,ntheta_constr
376 C ii = itheta_constr(i)
377 C thetabound(1,ii) = phi0(i)-drange(i)
378 C thetabound(2,ii) = phi0(i)+drange(i)
380 endif ! ntheta_constr.gt.0
381 endif! with_theta_constr
385 print *,'NNT=',NNT,' NCT=',NCT
386 if (itype(1).eq.ntyp1) nnt=2
387 if (itype(nres).eq.ntyp1) nct=nct-1
388 if (nstart.lt.nnt) nstart=nnt
389 if (nend.gt.nct .or. nend.eq.0) nend=nct
390 write (iout,*) "nstart",nstart," nend",nend
393 c read(inp,'(a)') pdbfile
394 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
395 c open(ipdbin,file=pdbfile,status='old',err=33)
397 c 33 write (iout,'(a)') 'Error opening PDB file.'
400 c print *,'Begin reading pdb data'
402 c print *,'Finished reading pdb data'
403 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
405 c itype_pdb(i)=itype(i)
408 c write (iout,'(a,i3)') 'nsup=',nsup
410 c if (nsup.le.(nct-nnt+1)) then
411 c do i=0,nct-nnt+1-nsup
412 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
418 c & 'Error - sequences to be superposed do not match.'
421 c do i=0,nsup-(nct-nnt+1)
422 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
424 c nstart_sup=nstart_sup+i
430 c & 'Error - sequences to be superposed do not match.'
433 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
434 c & ' nstart_seq=',nstart_seq
438 write (iout,*) "molread: REFSTR",refstr
440 if (.not.pdbref) then
441 call read_angles(inp,*38)
443 38 write (iout,'(a)') 'Error reading reference structure.'
445 call mp_stopall(Error_Msg)
447 stop 'Error reading reference structure'
460 call contact(.true.,ncont_ref,icont_ref)
463 C write (iout,'(/a,i3,a)')
464 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
465 write (iout,'(20i4)') (iss(i),i=1,ns)
467 write(iout,*)"Running with dynamic disulfide-bond formation"
469 write (iout,'(/a/)') 'Pre-formed links are:'
475 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
476 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
482 if (ns.gt.0.and.dyn_ss) then
486 forcon(i-nss)=forcon(i)
493 dyn_ss_mask(iss(i))=.true.
496 c Read distance restraints
497 if (constr_dist.gt.0) then
498 call read_dist_constr
503 c-----------------------------------------------------------------------------
504 logical function seq_comp(itypea,itypeb,length)
506 integer length,itypea(length),itypeb(length)
509 if (itypea(i).ne.itypeb(i)) then
517 c-----------------------------------------------------------------------------
518 subroutine read_bridge
519 C Read information about disulfide bridges.
522 include 'COMMON.IOUNITS'
525 include 'COMMON.INTERACT'
526 include 'COMMON.LOCAL'
527 include 'COMMON.NAMES'
528 include 'COMMON.CHAIN'
529 include 'COMMON.FFIELD'
530 include 'COMMON.SBRIDGE'
531 include 'COMMON.HEADER'
532 include 'COMMON.CONTROL'
533 include 'COMMON.TIME1'
535 include 'COMMON.INFO'
538 C Read bridging residues.
539 read (inp,*) ns,(iss(i),i=1,ns)
541 C Check whether the specified bridging residues are cystines.
543 if (itype(iss(i)).ne.1) then
544 write (iout,'(2a,i3,a)')
545 & 'Do you REALLY think that the residue ',
546 & restyp(itype(iss(i))),i,
547 & ' can form a disulfide bridge?!!!'
548 write (*,'(2a,i3,a)')
549 & 'Do you REALLY think that the residue ',
550 & restyp(itype(iss(i))),i,
551 & ' can form a disulfide bridge?!!!'
553 call mp_stopall(error_msg)
559 C Read preformed bridges.
561 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
564 C Check if the residues involved in bridges are in the specified list of
568 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
569 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
570 write (iout,'(a,i3,a)') 'Disulfide pair',i,
571 & ' contains residues present in other pairs.'
572 write (*,'(a,i3,a)') 'Disulfide pair',i,
573 & ' contains residues present in other pairs.'
575 call mp_stopall(error_msg)
582 if (ihpb(i).eq.iss(j)) goto 10
584 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
587 if (jhpb(i).eq.iss(j)) goto 20
589 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
602 c----------------------------------------------------------------------------
603 subroutine read_angles(kanal,*)
608 include 'COMMON.CHAIN'
609 include 'COMMON.IOUNITS'
611 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
612 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
613 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
614 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
616 theta(i)=deg2rad*theta(i)
617 phi(i)=deg2rad*phi(i)
618 alph(i)=deg2rad*alph(i)
619 omeg(i)=deg2rad*omeg(i)
624 c----------------------------------------------------------------------------
625 subroutine reada(rekord,lancuch,wartosc,default)
627 character*(*) rekord,lancuch
628 double precision wartosc,default
631 iread=index(rekord,lancuch)
636 iread=iread+ilen(lancuch)+1
637 read (rekord(iread:),*) wartosc
640 c----------------------------------------------------------------------------
641 subroutine multreada(rekord,lancuch,tablica,dim,default)
644 double precision tablica(dim),default
645 character*(*) rekord,lancuch
651 iread=index(rekord,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)
657 c----------------------------------------------------------------------------
658 subroutine readi(rekord,lancuch,wartosc,default)
660 character*(*) rekord,lancuch
661 integer wartosc,default
664 iread=index(rekord,lancuch)
669 iread=iread+ilen(lancuch)+1
670 read (rekord(iread:),*) wartosc
673 C----------------------------------------------------------------------
674 subroutine multreadi(rekord,lancuch,tablica,dim,default)
677 integer tablica(dim),default
678 character*(*) rekord,lancuch
685 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
686 if (iread.eq.0) return
687 iread=iread+ilen(lancuch)+1
688 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
692 c----------------------------------------------------------------------------
693 subroutine card_concat(card)
695 include 'COMMON.IOUNITS'
697 character*80 karta,ucase
699 read (inp,'(a)') karta
702 do while (karta(80:80).eq.'&')
703 card=card(:ilen(card)+1)//karta(:79)
704 read (inp,'(a)') karta
707 card=card(:ilen(card)+1)//karta
710 c----------------------------------------------------------------------------
719 include 'COMMON.IOUNITS'
720 include 'COMMON.CONTROL'
721 integer lenpre,lenpot,ilen
723 character*16 cformat,cprint
725 integer lenint,lenout
726 call getenv('INPUT',prefix)
727 call getenv('OUTPUT',prefout)
728 call getenv('INTIN',prefintin)
729 call getenv('COORD',cformat)
730 call getenv('PRINTCOOR',cprint)
731 call getenv('SCRATCHDIR',scratchdir)
734 if (index(ucase(cformat),'CX').gt.0) then
741 lenint=ilen(prefintin)
742 C Get the names and open the input files
743 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
745 write (liczba,'(bz,i3.3)') me
746 outname=prefout(:lenout)//'_clust.out_'//liczba
748 outname=prefout(:lenout)//'_clust.out'
751 intinname=prefintin(:lenint)//'.bx'
752 else if (from_cx) then
753 intinname=prefintin(:lenint)//'.cx'
755 intinname=prefintin(:lenint)//'.int'
757 rmsname=prefintin(:lenint)//'.rms'
758 open (jplot,file=prefout(:ilen(prefout))//'.tex',
760 open (jrms,file=rmsname,status='unknown')
761 open(iout,file=outname,status='unknown')
762 C Get parameter filenames and open the parameter files.
763 call getenv('BONDPAR',bondname)
764 open (ibond,file=bondname,status='old')
765 call getenv('THETPAR',thetname)
766 open (ithep,file=thetname,status='old')
767 call getenv('ROTPAR',rotname)
768 open (irotam,file=rotname,status='old')
769 call getenv('TORPAR',torname)
770 open (itorp,file=torname,status='old')
771 call getenv('TORDPAR',tordname)
772 open (itordp,file=tordname,status='old')
773 call getenv('FOURIER',fouriername)
774 open (ifourier,file=fouriername,status='old')
775 call getenv('ELEPAR',elename)
776 open (ielep,file=elename,status='old')
777 call getenv('SIDEPAR',sidename)
778 open (isidep,file=sidename,status='old')
779 call getenv('SIDEP',sidepname)
780 open (isidep1,file=sidepname,status="old")
781 call getenv('SCCORPAR',sccorname)
782 open (isccor,file=sccorname,status="old")
785 C 8/9/01 In the newest version SCp interaction constants are read from a file
786 C Use -DOLDSCP to use hard-coded constants instead.
788 call getenv('SCPPAR',scpname)
789 open (iscpp,file=scpname,status='old')
793 subroutine read_dist_constr
794 implicit real*8 (a-h,o-z)
799 include 'COMMON.CONTROL'
800 include 'COMMON.CHAIN'
801 include 'COMMON.IOUNITS'
802 include 'COMMON.SBRIDGE'
803 integer ifrag_(2,100),ipair_(2,100)
804 double precision wfrag_(100),wpair_(100)
805 character*500 controlcard
806 logical lprn /.true./
807 write (iout,*) "Calling read_dist_constr"
808 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
810 write(iout,*) "TU sie wywalam?"
811 call card_concat(controlcard)
812 write (iout,*) controlcard
814 call readi(controlcard,"NFRAG",nfrag_,0)
815 call readi(controlcard,"NPAIR",npair_,0)
816 call readi(controlcard,"NDIST",ndist_,0)
817 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
818 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
819 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
820 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
821 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
822 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
823 write (iout,*) "IFRAG"
825 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
827 write (iout,*) "IPAIR"
829 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
833 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
834 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
835 & ifrag_(2,i)=nstart_sup+nsup-1
836 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
838 if (wfrag_(i).gt.0.0d0) then
839 do j=ifrag_(1,i),ifrag_(2,i)-1
841 write (iout,*) "j",j," k",k
843 if (constr_dist.eq.1) then
848 forcon(nhpb)=wfrag_(i)
849 else if (constr_dist.eq.2) then
850 if (ddjk.le.dist_cut) then
855 forcon(nhpb)=wfrag_(i)
862 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
865 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
866 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
872 if (wpair_(i).gt.0.0d0) then
880 do j=ifrag_(1,ii),ifrag_(2,ii)
881 do k=ifrag_(1,jj),ifrag_(2,jj)
885 forcon(nhpb)=wpair_(i)
887 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
888 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
894 if (constr_dist.eq.11) then
895 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
896 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
897 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
898 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
899 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
901 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
903 if (forcon(nhpb+1).gt.0.0d0) then
905 if (ibecarb(i).gt.0) then
909 if (dhpb(nhpb).eq.0.0d0)
910 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
911 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
912 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
913 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)