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)
40 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
41 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
42 if (lipthick.gt.0.0d0) then
43 bordliptop=(boxzsize+lipthick)/2.0
44 bordlipbot=bordliptop-lipthick
46 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
47 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
48 buflipbot=bordlipbot+lipbufthick
49 bufliptop=bordliptop-lipbufthick
50 if ((lipbufthick*2.0d0).gt.lipthick)
51 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
53 write(iout,*) "bordliptop=",bordliptop
54 write(iout,*) "bordlipbot=",bordlipbot
55 write(iout,*) "bufliptop=",bufliptop
56 write(iout,*) "buflipbot=",buflipbot
58 call readi(controlcard,'SHIELD',shield_mode,0)
59 write (iout,*) "SHIELD MODE",shield_mode
60 call readi(controlcard,'TUBEMOD',tubelog,0)
61 write (iout,*) TUBElog,"TUBEMODE"
62 if (shield_mode.gt.0) then
64 C VSolvSphere the volume of solving sphere
66 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
67 C there will be no distinction between proline peptide group and normal peptide
68 C group in case of shielding parameters
69 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
70 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
71 write (iout,*) VSolvSphere,VSolvSphere_div
72 C long axis of side chain
74 long_r_sidechain(i)=vbldsc0(1,i)
75 short_r_sidechain(i)=sigma0(i)
79 if (TUBElog.gt.0) then
80 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
81 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
82 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
83 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
84 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
85 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
86 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
87 buftubebot=bordtubebot+tubebufthick
88 buftubetop=bordtubetop-tubebufthick
90 call readi(controlcard,'PDBOUT',outpdb,0)
91 call readi(controlcard,'MOL2OUT',outmol2,0)
92 refstr=(index(controlcard,'REFSTR').gt.0)
93 write (iout,*) "REFSTR",refstr
94 pdbref=(index(controlcard,'PDBREF').gt.0)
95 iscode=index(controlcard,'ONE_LETTER')
96 tree=(index(controlcard,'MAKE_TREE').gt.0)
97 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 write (iout,*) "with_dihed_constr ",with_dihed_constr,
100 & " CONSTR_DIST",constr_dist
101 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
102 write (iout,*) "with_theta_constr ",with_theta_constr
104 min_var=(index(controlcard,'MINVAR').gt.0)
105 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
106 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
107 call readi(controlcard,'NCUT',ncut,1)
108 call readi(controlcard,'SYM',symetr,1)
109 write (iout,*) 'sym', symetr
110 call readi(controlcard,'NSTART',nstart,0)
111 call readi(controlcard,'NEND',nend,0)
112 call reada(controlcard,'ECUT',ecut,10.0d0)
113 call reada(controlcard,'PROB',prob_limit,0.99d0)
114 write (iout,*) "Probability limit",prob_limit
115 lgrp=(index(controlcard,'LGRP').gt.0)
116 caonly=(index(controlcard,'CA_ONLY').gt.0)
117 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
118 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
119 call readi(controlcard,'IOPT',iopt,2)
120 lside = index(controlcard,"SIDE").gt.0
121 efree = index(controlcard,"EFREE").gt.0
122 call readi(controlcard,'NTEMP',nT,1)
123 write (iout,*) "nT",nT
124 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
125 write (iout,*) "nT",nT
126 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
128 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
130 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
131 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
132 lprint_int=index(controlcard,"PRINT_INT") .gt.0
136 c--------------------------------------------------------------------------
139 C Read molecular data.
143 include 'COMMON.IOUNITS'
146 include 'COMMON.INTERACT'
147 include 'COMMON.LOCAL'
148 include 'COMMON.NAMES'
149 include 'COMMON.CHAIN'
150 include 'COMMON.FFIELD'
151 include 'COMMON.SBRIDGE'
152 include 'COMMON.HEADER'
153 include 'COMMON.CONTROL'
154 include 'COMMON.CONTACTS'
155 include 'COMMON.TIME1'
156 include 'COMMON.TORCNSTR'
157 include 'COMMON.SHIELD'
159 include 'COMMON.INFO'
161 character*4 sequence(maxres)
162 character*800 weightcard
164 double precision x(maxvar)
165 integer itype_pdb(maxres)
167 integer i,j,kkk,i1,i2,it1,it2
171 C Read weights of the subsequent energy terms.
172 call card_concat(weightcard)
173 write(iout,*) weightcard
174 C call reada(weightcard,'WSC',wsc,1.0d0)
176 call reada(weightcard,'WLONG',wsc,wsc)
177 call reada(weightcard,'WSCP',wscp,1.0d0)
178 call reada(weightcard,'WELEC',welec,1.0D0)
179 call reada(weightcard,'WVDWPP',wvdwpp,welec)
180 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
181 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
182 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
183 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
184 call reada(weightcard,'WTURN3',wturn3,1.0D0)
185 call reada(weightcard,'WTURN4',wturn4,1.0D0)
186 call reada(weightcard,'WTURN6',wturn6,1.0D0)
187 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
188 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
189 call reada(weightcard,'WBOND',wbond,1.0D0)
190 call reada(weightcard,'WTOR',wtor,1.0D0)
191 call reada(weightcard,'WTORD',wtor_d,1.0D0)
192 call reada(weightcard,'WANG',wang,1.0D0)
193 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
194 call reada(weightcard,'SCAL14',scal14,0.4D0)
195 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
196 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
197 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
198 if (index(weightcard,'SOFT').gt.0) ipot=6
199 call reada(weightcard,"D0CM",d0cm,3.78d0)
200 call reada(weightcard,"AKCM",akcm,15.1d0)
201 call reada(weightcard,"AKTH",akth,11.0d0)
202 call reada(weightcard,"AKCT",akct,12.0d0)
203 call reada(weightcard,"V1SS",v1ss,-1.08d0)
204 call reada(weightcard,"V2SS",v2ss,7.61d0)
205 call reada(weightcard,"V3SS",v3ss,13.7d0)
206 call reada(weightcard,"EBR",ebr,-5.50D0)
207 call reada(weightcard,'WSHIELD',wshield,1.0d0)
208 write(iout,*) 'WSHIELD',wshield
209 call reada(weightcard,'WTUBE',wtube,0.0d0)
210 write(iout,*) 'WTUBE',wtube
211 call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
212 call reada(weightcard,'WLT',wliptran,0.0D0)
213 call reada(weightcard,"ATRISS",atriss,0.301D0)
214 call reada(weightcard,"BTRISS",btriss,0.021D0)
215 call reada(weightcard,"CTRISS",ctriss,1.001D0)
216 call reada(weightcard,"DTRISS",dtriss,1.001D0)
217 write (iout,*) "ATRISS=", atriss
218 write (iout,*) "BTRISS=", btriss
219 write (iout,*) "CTRISS=", ctriss
220 write (iout,*) "DTRISS=", dtriss
221 if (shield_mode.gt.0) then
223 write (iout,*) "WARNING:liscale not used in shield mode"
225 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
227 dyn_ss_mask(i)=.false.
231 dyn_ssbond_ij(i,j)=1.0d300
234 call reada(weightcard,"HT",Ht,0.0D0)
236 ss_depth=ebr/wsc-0.25*eps(1,1)
237 Ht=Ht/wsc-0.25*eps(1,1)
238 akcm=akcm*wstrain/wsc
239 akth=akth*wstrain/wsc
240 akct=akct*wstrain/wsc
241 v1ss=v1ss*wstrain/wsc
242 v2ss=v2ss*wstrain/wsc
243 v3ss=v3ss*wstrain/wsc
245 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
247 write (iout,'(/a)') "Disulfide bridge parameters:"
248 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
249 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
250 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
251 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
254 C 12/1/95 Added weight for the multi-body term WCORR
255 call reada(weightcard,'WCORRH',wcorr,1.0D0)
256 if (wcorr4.gt.0.0d0) wcorr=wcorr4
276 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
277 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
278 & wturn4,wturn6,wsccor
279 10 format (/'Energy-term weights (unscaled):'//
280 & 'WSCC= ',f10.6,' (SC-SC)'/
281 & 'WSCP= ',f10.6,' (SC-p)'/
282 & 'WELEC= ',f10.6,' (p-p electr)'/
283 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
284 & 'WBOND= ',f10.6,' (stretching)'/
285 & 'WANG= ',f10.6,' (bending)'/
286 & 'WSCLOC= ',f10.6,' (SC local)'/
287 & 'WTOR= ',f10.6,' (torsional)'/
288 & 'WTORD= ',f10.6,' (double torsional)'/
289 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
290 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
291 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
292 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
293 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
294 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
295 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
296 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
297 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
299 if (wcorr4.gt.0.0d0) then
300 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
301 & 'between contact pairs of peptide groups'
302 write (iout,'(2(a,f5.3/))')
303 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
304 & 'Range of quenching the correlation terms:',2*delt_corr
305 else if (wcorr.gt.0.0d0) then
306 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
307 & 'between contact pairs of peptide groups'
309 write (iout,'(a,f8.3)')
310 & 'Scaling factor of 1,4 SC-p interactions:',scal14
311 write (iout,'(a,f8.3)')
312 & 'General scaling factor of SC-p interactions:',scalscp
313 r0_corr=cutoff_corr-delt_corr
315 aad(i,1)=scalscp*aad(i,1)
316 aad(i,2)=scalscp*aad(i,2)
317 bad(i,1)=scalscp*bad(i,1)
318 bad(i,2)=scalscp*bad(i,2)
322 print *,'indpdb=',indpdb,' pdbref=',pdbref
324 C Read sequence if not taken from the pdb file.
325 if (iscode.gt.0) then
326 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
328 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
330 C Convert sequence to numeric code
332 itype(i)=rescode(i,sequence(i),iscode)
335 print '(20i4)',(itype(i),i=1,nres)
339 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
341 if (itype(i).eq.ntyp1) then
345 else if (iabs(itype(i+1)).ne.20) then
347 else if (iabs(itype(i)).ne.20) then
354 write (iout,*) "ITEL"
356 write (iout,*) i,itype(i),itel(i)
359 print *,'Call Read_Bridge.'
361 C this fragment reads diheadral constrains
362 if (with_dihed_constr) then
364 read (inp,*) ndih_constr
365 if (ndih_constr.gt.0) then
367 C write (iout,*) 'FTORS',ftors
368 C ftors is the force constant for torsional quartic constrains
369 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
372 & 'There are',ndih_constr,' constraints on phi angles.'
374 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
378 phi0(i)=deg2rad*phi0(i)
379 drange(i)=deg2rad*drange(i)
381 endif ! endif ndif_constr.gt.0
382 endif ! with_dihed_constr
383 if (with_theta_constr) then
384 C with_theta_constr is keyword allowing for occurance of theta constrains
385 read (inp,*) ntheta_constr
386 C ntheta_constr is the number of theta constrains
387 if (ntheta_constr.gt.0) then
389 read (inp,*) (itheta_constr(i),theta_constr0(i),
390 & theta_drange(i),for_thet_constr(i),
392 C the above code reads from 1 to ntheta_constr
393 C itheta_constr(i) residue i for which is theta_constr
394 C theta_constr0 the global minimum value
395 C theta_drange is range for which there is no energy penalty
396 C for_thet_constr is the force constant for quartic energy penalty
398 C if(me.eq.king.or..not.out1file)then
400 & 'There are',ntheta_constr,' constraints on phi angles.'
402 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
408 theta_constr0(i)=deg2rad*theta_constr0(i)
409 theta_drange(i)=deg2rad*theta_drange(i)
411 C if(me.eq.king.or..not.out1file)
412 C & write (iout,*) 'FTORS',ftors
413 C do i=1,ntheta_constr
414 C ii = itheta_constr(i)
415 C thetabound(1,ii) = phi0(i)-drange(i)
416 C thetabound(2,ii) = phi0(i)+drange(i)
418 endif ! ntheta_constr.gt.0
419 endif! with_theta_constr
423 print *,'NNT=',NNT,' NCT=',NCT
424 if (itype(1).eq.ntyp1) nnt=2
425 if (itype(nres).eq.ntyp1) nct=nct-1
426 if (nstart.lt.nnt) nstart=nnt
427 if (nend.gt.nct .or. nend.eq.0) nend=nct
428 write (iout,*) "nstart",nstart," nend",nend
431 c read(inp,'(a)') pdbfile
432 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
433 c open(ipdbin,file=pdbfile,status='old',err=33)
435 c 33 write (iout,'(a)') 'Error opening PDB file.'
438 c print *,'Begin reading pdb data'
440 c print *,'Finished reading pdb data'
441 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
443 c itype_pdb(i)=itype(i)
446 c write (iout,'(a,i3)') 'nsup=',nsup
448 c if (nsup.le.(nct-nnt+1)) then
449 c do i=0,nct-nnt+1-nsup
450 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
456 c & 'Error - sequences to be superposed do not match.'
459 c do i=0,nsup-(nct-nnt+1)
460 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
462 c nstart_sup=nstart_sup+i
468 c & 'Error - sequences to be superposed do not match.'
471 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
472 c & ' nstart_seq=',nstart_seq
476 write (iout,*) "molread: REFSTR",refstr
478 if (.not.pdbref) then
479 call read_angles(inp,*38)
481 38 write (iout,'(a)') 'Error reading reference structure.'
483 call mp_stopall(Error_Msg)
485 stop 'Error reading reference structure'
498 call contact(.true.,ncont_ref,icont_ref)
501 C write (iout,'(/a,i3,a)')
502 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
503 write (iout,'(20i4)') (iss(i),i=1,ns)
505 write(iout,*)"Running with dynamic disulfide-bond formation"
507 write (iout,'(/a/)') 'Pre-formed links are:'
513 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
514 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
520 if (ns.gt.0.and.dyn_ss) then
524 forcon(i-nss)=forcon(i)
531 dyn_ss_mask(iss(i))=.true.
534 c Read distance restraints
535 if (constr_dist.gt.0) then
536 call read_dist_constr
541 c-----------------------------------------------------------------------------
542 logical function seq_comp(itypea,itypeb,length)
544 integer length,itypea(length),itypeb(length)
547 if (itypea(i).ne.itypeb(i)) then
555 c-----------------------------------------------------------------------------
556 subroutine read_bridge
557 C Read information about disulfide bridges.
560 include 'COMMON.IOUNITS'
563 include 'COMMON.INTERACT'
564 include 'COMMON.LOCAL'
565 include 'COMMON.NAMES'
566 include 'COMMON.CHAIN'
567 include 'COMMON.FFIELD'
568 include 'COMMON.SBRIDGE'
569 include 'COMMON.HEADER'
570 include 'COMMON.CONTROL'
571 include 'COMMON.TIME1'
573 include 'COMMON.INFO'
576 C Read bridging residues.
577 read (inp,*) ns,(iss(i),i=1,ns)
579 C Check whether the specified bridging residues are cystines.
581 if (itype(iss(i)).ne.1) then
582 write (iout,'(2a,i3,a)')
583 & 'Do you REALLY think that the residue ',
584 & restyp(itype(iss(i))),i,
585 & ' can form a disulfide bridge?!!!'
586 write (*,'(2a,i3,a)')
587 & 'Do you REALLY think that the residue ',
588 & restyp(itype(iss(i))),i,
589 & ' can form a disulfide bridge?!!!'
591 call mp_stopall(error_msg)
597 C Read preformed bridges.
599 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
602 C Check if the residues involved in bridges are in the specified list of
606 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
607 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
608 write (iout,'(a,i3,a)') 'Disulfide pair',i,
609 & ' contains residues present in other pairs.'
610 write (*,'(a,i3,a)') 'Disulfide pair',i,
611 & ' contains residues present in other pairs.'
613 call mp_stopall(error_msg)
620 if (ihpb(i).eq.iss(j)) goto 10
622 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
625 if (jhpb(i).eq.iss(j)) goto 20
627 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
640 c----------------------------------------------------------------------------
641 subroutine read_angles(kanal,*)
646 include 'COMMON.CHAIN'
647 include 'COMMON.IOUNITS'
649 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
650 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
651 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
652 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
654 theta(i)=deg2rad*theta(i)
655 phi(i)=deg2rad*phi(i)
656 alph(i)=deg2rad*alph(i)
657 omeg(i)=deg2rad*omeg(i)
662 c----------------------------------------------------------------------------
663 subroutine reada(rekord,lancuch,wartosc,default)
665 character*(*) rekord,lancuch
666 double precision wartosc,default
669 iread=index(rekord,lancuch)
674 iread=iread+ilen(lancuch)+1
675 read (rekord(iread:),*) wartosc
678 c----------------------------------------------------------------------------
679 subroutine multreada(rekord,lancuch,tablica,dim,default)
682 double precision tablica(dim),default
683 character*(*) rekord,lancuch
689 iread=index(rekord,lancuch)
690 if (iread.eq.0) return
691 iread=iread+ilen(lancuch)+1
692 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
695 c----------------------------------------------------------------------------
696 subroutine readi(rekord,lancuch,wartosc,default)
698 character*(*) rekord,lancuch
699 integer wartosc,default
702 iread=index(rekord,lancuch)
707 iread=iread+ilen(lancuch)+1
708 read (rekord(iread:),*) wartosc
711 C----------------------------------------------------------------------
712 subroutine multreadi(rekord,lancuch,tablica,dim,default)
715 integer tablica(dim),default
716 character*(*) rekord,lancuch
723 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
724 if (iread.eq.0) return
725 iread=iread+ilen(lancuch)+1
726 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
730 c----------------------------------------------------------------------------
731 subroutine card_concat(card)
733 include 'COMMON.IOUNITS'
735 character*80 karta,ucase
737 read (inp,'(a)') karta
740 do while (karta(80:80).eq.'&')
741 card=card(:ilen(card)+1)//karta(:79)
742 read (inp,'(a)') karta
745 card=card(:ilen(card)+1)//karta
748 c----------------------------------------------------------------------------
757 include 'COMMON.IOUNITS'
758 include 'COMMON.CONTROL'
759 integer lenpre,lenpot,ilen
761 character*16 cformat,cprint
763 integer lenint,lenout
764 call getenv('INPUT',prefix)
765 call getenv('OUTPUT',prefout)
766 call getenv('INTIN',prefintin)
767 call getenv('COORD',cformat)
768 call getenv('PRINTCOOR',cprint)
769 call getenv('SCRATCHDIR',scratchdir)
772 if (index(ucase(cformat),'CX').gt.0) then
779 lenint=ilen(prefintin)
780 C Get the names and open the input files
781 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
783 write (liczba,'(bz,i3.3)') me
784 outname=prefout(:lenout)//'_clust.out_'//liczba
786 outname=prefout(:lenout)//'_clust.out'
789 intinname=prefintin(:lenint)//'.bx'
790 else if (from_cx) then
791 intinname=prefintin(:lenint)//'.cx'
793 intinname=prefintin(:lenint)//'.int'
795 rmsname=prefintin(:lenint)//'.rms'
796 open (jplot,file=prefout(:ilen(prefout))//'.tex',
798 open (jrms,file=rmsname,status='unknown')
799 open(iout,file=outname,status='unknown')
800 C Get parameter filenames and open the parameter files.
801 call getenv('BONDPAR',bondname)
802 open (ibond,file=bondname,status='old')
803 call getenv('THETPAR',thetname)
804 open (ithep,file=thetname,status='old')
805 call getenv('ROTPAR',rotname)
806 open (irotam,file=rotname,status='old')
807 call getenv('TORPAR',torname)
808 open (itorp,file=torname,status='old')
809 call getenv('TORDPAR',tordname)
810 open (itordp,file=tordname,status='old')
811 call getenv('FOURIER',fouriername)
812 open (ifourier,file=fouriername,status='old')
813 call getenv('ELEPAR',elename)
814 open (ielep,file=elename,status='old')
815 call getenv('SIDEPAR',sidename)
816 open (isidep,file=sidename,status='old')
817 call getenv('SIDEP',sidepname)
818 open (isidep1,file=sidepname,status="old")
819 call getenv('SCCORPAR',sccorname)
820 open (isccor,file=sccorname,status="old")
821 call getenv('LIPTRANPAR',liptranname)
822 open (iliptranpar,file=liptranname,status='old')
823 call getenv('TUBEPAR',tubename)
824 open (itube,file=tubename,status='old')
828 C 8/9/01 In the newest version SCp interaction constants are read from a file
829 C Use -DOLDSCP to use hard-coded constants instead.
831 call getenv('SCPPAR',scpname)
832 open (iscpp,file=scpname,status='old')
836 subroutine read_dist_constr
837 implicit real*8 (a-h,o-z)
842 include 'COMMON.CONTROL'
843 include 'COMMON.CHAIN'
844 include 'COMMON.IOUNITS'
845 include 'COMMON.SBRIDGE'
846 integer ifrag_(2,100),ipair_(2,100)
847 double precision wfrag_(100),wpair_(100)
848 character*500 controlcard
849 logical lprn /.true./
850 write (iout,*) "Calling read_dist_constr"
851 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
853 write(iout,*) "TU sie wywalam?"
854 call card_concat(controlcard)
855 write (iout,*) controlcard
857 call readi(controlcard,"NFRAG",nfrag_,0)
858 call readi(controlcard,"NPAIR",npair_,0)
859 call readi(controlcard,"NDIST",ndist_,0)
860 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
861 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
862 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
863 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
864 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
865 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
866 write (iout,*) "IFRAG"
868 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
870 write (iout,*) "IPAIR"
872 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
876 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
877 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
878 & ifrag_(2,i)=nstart_sup+nsup-1
879 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
881 if (wfrag_(i).gt.0.0d0) then
882 do j=ifrag_(1,i),ifrag_(2,i)-1
884 write (iout,*) "j",j," k",k
886 if (constr_dist.eq.1) then
891 forcon(nhpb)=wfrag_(i)
892 else if (constr_dist.eq.2) then
893 if (ddjk.le.dist_cut) then
898 forcon(nhpb)=wfrag_(i)
905 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
908 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
909 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
915 if (wpair_(i).gt.0.0d0) then
923 do j=ifrag_(1,ii),ifrag_(2,ii)
924 do k=ifrag_(1,jj),ifrag_(2,jj)
928 forcon(nhpb)=wpair_(i)
930 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
931 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
937 if (constr_dist.eq.11) then
938 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
939 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
940 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
941 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
942 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
944 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
946 if (forcon(nhpb+1).gt.0.0d0) then
948 if (ibecarb(i).gt.0) then
952 if (dhpb(nhpb).eq.0.0d0)
953 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
954 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
955 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
956 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)