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 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
30 call readi(controlcard,'TORMODE',tor_mode,0)
31 call readi(controlcard,'NRES',nres,0)
32 call readi(controlcard,'RESCALE',rescale_mode,2)
33 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
34 write (iout,*) "DISTCHAINMAX",distchainmax
35 C Reading the dimensions of box in x,y,z coordinates
36 call reada(controlcard,'BOXX',boxxsize,100.0d0)
37 call reada(controlcard,'BOXY',boxysize,100.0d0)
38 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
39 c Cutoff range for interactions
40 call reada(controlcard,"R_CUT",r_cut,15.0d0)
41 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
42 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
43 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
44 if (lipthick.gt.0.0d0) then
45 bordliptop=(boxzsize+lipthick)/2.0
46 bordlipbot=bordliptop-lipthick
48 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
49 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
50 buflipbot=bordlipbot+lipbufthick
51 bufliptop=bordliptop-lipbufthick
52 if ((lipbufthick*2.0d0).gt.lipthick)
53 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
55 write(iout,*) "bordliptop=",bordliptop
56 write(iout,*) "bordlipbot=",bordlipbot
57 write(iout,*) "bufliptop=",bufliptop
58 write(iout,*) "buflipbot=",buflipbot
60 call readi(controlcard,'SHIELD',shield_mode,0)
61 write (iout,*) "SHIELD MODE",shield_mode
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 call readi(controlcard,'PDBOUT',outpdb,0)
80 call readi(controlcard,'MOL2OUT',outmol2,0)
81 refstr=(index(controlcard,'REFSTR').gt.0)
82 pdbref=(index(controlcard,'PDBREF').gt.0)
83 refstr = refstr .or. pdbref
84 write (iout,*) "REFSTR",refstr," PDBREF",pdbref
85 iscode=index(controlcard,'ONE_LETTER')
86 tree=(index(controlcard,'MAKE_TREE').gt.0)
87 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
88 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
89 write (iout,*) "with_dihed_constr ",with_dihed_constr,
90 & " CONSTR_DIST",constr_dist
91 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
92 write (iout,*) "with_theta_constr ",with_theta_constr
94 min_var=(index(controlcard,'MINVAR').gt.0)
95 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
96 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
97 call readi(controlcard,'NCUT',ncut,0)
99 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
102 call readi(controlcard,'NCLUST',nclust,5)
104 call readi(controlcard,'SYM',symetr,1)
105 write (iout,*) 'sym', symetr
106 call readi(controlcard,'NSTART',nstart,0)
107 call readi(controlcard,'NEND',nend,0)
108 call reada(controlcard,'ECUT',ecut,10.0d0)
109 call reada(controlcard,'PROB',prob_limit,0.99d0)
110 write (iout,*) "Probability limit",prob_limit
111 lgrp=(index(controlcard,'LGRP').gt.0)
112 caonly=(index(controlcard,'CA_ONLY').gt.0)
113 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
114 call readi(controlcard,'IOPT',iopt,2)
115 lside = index(controlcard,"SIDE").gt.0
116 lefree = index(controlcard,"EFREE").gt.0
117 call readi(controlcard,'NTEMP',nT,1)
118 write (iout,*) "nT",nT
119 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
120 write (iout,*) "nT",nT
121 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
123 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
125 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
126 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
127 lprint_int=index(controlcard,"PRINT_INT") .gt.0
131 c--------------------------------------------------------------------------
134 C Read molecular data.
138 include 'COMMON.IOUNITS'
141 include 'COMMON.INTERACT'
142 include 'COMMON.LOCAL'
143 include 'COMMON.NAMES'
144 include 'COMMON.CHAIN'
145 include 'COMMON.FFIELD'
146 include 'COMMON.SBRIDGE'
147 include 'COMMON.HEADER'
148 include 'COMMON.CONTROL'
149 include 'COMMON.CONTACTS'
150 include 'COMMON.TIME1'
151 include 'COMMON.TORCNSTR'
152 include 'COMMON.SHIELD'
154 include 'COMMON.INFO'
156 character*4 sequence(maxres)
157 character*800 weightcard,controlcard
159 double precision x(maxvar)
160 double precision phihel,phibet,sigmahel,sigmabet,sumv,
162 integer itype_pdb(maxres)
164 integer i,j,kkk,i1,i2,it1,it2
168 C Read weights of the subsequent energy terms.
169 call card_concat(weightcard)
170 call reada(weightcard,'WSC',wsc,1.0d0)
171 call reada(weightcard,'WLONG',wsc,wsc)
172 call reada(weightcard,'WSCP',wscp,1.0d0)
173 call reada(weightcard,'WELEC',welec,1.0D0)
174 call reada(weightcard,'WVDWPP',wvdwpp,welec)
175 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
176 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
177 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
178 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
179 call reada(weightcard,'WTURN3',wturn3,1.0D0)
180 call reada(weightcard,'WTURN4',wturn4,1.0D0)
181 call reada(weightcard,'WTURN6',wturn6,1.0D0)
182 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
183 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
184 call reada(weightcard,'WBOND',wbond,1.0D0)
185 call reada(weightcard,'WTOR',wtor,1.0D0)
186 call reada(weightcard,'WTORD',wtor_d,1.0D0)
187 call reada(weightcard,'WANG',wang,1.0D0)
188 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
189 call reada(weightcard,'SCAL14',scal14,0.4D0)
190 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
191 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
192 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
193 if (index(weightcard,'SOFT').gt.0) ipot=6
194 call reada(weightcard,"D0CM",d0cm,3.78d0)
195 call reada(weightcard,"AKCM",akcm,15.1d0)
196 call reada(weightcard,"AKTH",akth,11.0d0)
197 call reada(weightcard,"AKCT",akct,12.0d0)
198 call reada(weightcard,"V1SS",v1ss,-1.08d0)
199 call reada(weightcard,"V2SS",v2ss,7.61d0)
200 call reada(weightcard,"V3SS",v3ss,13.7d0)
201 call reada(weightcard,"EBR",ebr,-5.50D0)
202 call reada(weightcard,'WSHIELD',wshield,1.0d0)
203 write(iout,*) 'WSHIELD',wshield
204 call reada(weightcard,'WLT',wliptran,0.0D0)
205 call reada(weightcard,"ATRISS",atriss,0.301D0)
206 call reada(weightcard,"BTRISS",btriss,0.021D0)
207 call reada(weightcard,"CTRISS",ctriss,1.001D0)
208 call reada(weightcard,"DTRISS",dtriss,1.001D0)
209 write (iout,*) "ATRISS=", atriss
210 write (iout,*) "BTRISS=", btriss
211 write (iout,*) "CTRISS=", ctriss
212 write (iout,*) "DTRISS=", dtriss
213 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
215 dyn_ss_mask(i)=.false.
219 dyn_ssbond_ij(i,j)=1.0d300
222 call reada(weightcard,"HT",Ht,0.0D0)
224 ss_depth=ebr/wsc-0.25*eps(1,1)
225 Ht=Ht/wsc-0.25*eps(1,1)
226 akcm=akcm*wstrain/wsc
227 akth=akth*wstrain/wsc
228 akct=akct*wstrain/wsc
229 v1ss=v1ss*wstrain/wsc
230 v2ss=v2ss*wstrain/wsc
231 v3ss=v3ss*wstrain/wsc
233 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
235 write (iout,'(/a)') "Disulfide bridge parameters:"
236 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
237 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
238 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
239 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
242 C 12/1/95 Added weight for the multi-body term WCORR
243 call reada(weightcard,'WCORRH',wcorr,1.0D0)
244 if (wcorr4.gt.0.0d0) wcorr=wcorr4
264 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
265 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
266 & wturn4,wturn6,wsccor
267 10 format (/'Energy-term weights (unscaled):'//
268 & 'WSCC= ',f10.6,' (SC-SC)'/
269 & 'WSCP= ',f10.6,' (SC-p)'/
270 & 'WELEC= ',f10.6,' (p-p electr)'/
271 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
272 & 'WBOND= ',f10.6,' (stretching)'/
273 & 'WANG= ',f10.6,' (bending)'/
274 & 'WSCLOC= ',f10.6,' (SC local)'/
275 & 'WTOR= ',f10.6,' (torsional)'/
276 & 'WTORD= ',f10.6,' (double torsional)'/
277 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
278 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
279 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
280 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
281 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
282 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
283 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
284 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
285 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
287 if (wcorr4.gt.0.0d0) then
288 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
289 & 'between contact pairs of peptide groups'
290 write (iout,'(2(a,f5.3/))')
291 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
292 & 'Range of quenching the correlation terms:',2*delt_corr
293 else if (wcorr.gt.0.0d0) then
294 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
295 & 'between contact pairs of peptide groups'
297 write (iout,'(a,f8.3)')
298 & 'Scaling factor of 1,4 SC-p interactions:',scal14
299 write (iout,'(a,f8.3)')
300 & 'General scaling factor of SC-p interactions:',scalscp
301 r0_corr=cutoff_corr-delt_corr
303 aad(i,1)=scalscp*aad(i,1)
304 aad(i,2)=scalscp*aad(i,2)
305 bad(i,1)=scalscp*bad(i,1)
306 bad(i,2)=scalscp*bad(i,2)
310 c print *,'indpdb=',indpdb,' pdbref=',pdbref
312 C Read sequence if not taken from the pdb file.
313 if (iscode.gt.0) then
314 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
316 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
318 C Convert sequence to numeric code
320 itype(i)=rescode(i,sequence(i),iscode)
323 c print '(20i4)',(itype(i),i=1,nres)
327 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
329 if (itype(i).eq.ntyp1) then
333 else if (iabs(itype(i+1)).ne.20) then
335 else if (iabs(itype(i)).ne.20) then
342 write (iout,*) "ITEL"
344 write (iout,*) i,itype(i),itel(i)
347 c print *,'Call Read_Bridge.'
349 C this fragment reads diheadral constrains
352 c print *,'NNT=',NNT,' NCT=',NCT
353 if (itype(1).eq.ntyp1) nnt=2
354 if (itype(nres).eq.ntyp1) nct=nct-1
355 if (nstart.lt.nnt) nstart=nnt
356 if (nend.gt.nct .or. nend.eq.0) nend=nct
357 write (iout,*) "nstart",nstart," nend",nend
359 if (with_dihed_constr) then
361 read (inp,*) ndih_constr
362 if (ndih_constr.gt.0) then
365 C write (iout,*) 'FTORS',ftors
366 C ftors is the force constant for torsional quartic constrains
367 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
370 & 'There are',ndih_constr,' constraints on phi angles.'
372 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
376 phi0(i)=deg2rad*phi0(i)
377 drange(i)=deg2rad*drange(i)
379 else if (ndih_constr.lt.0) then
381 call card_concat(controlcard)
382 call reada(controlcard,"PHIHEL",phihel,50.0D0)
383 call reada(controlcard,"PHIBET",phibet,180.0D0)
384 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
385 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
386 call reada(controlcard,"WDIHC",wdihc,0.591d0)
387 write (iout,*) "Weight of the dihedral restraint term",wdihc
388 read(inp,'(9x,3f7.3)')
389 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
390 write (iout,*) "The secprob array"
392 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
396 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
397 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
398 ndih_constr=ndih_constr+1
399 idih_constr(ndih_constr)=i
402 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
403 sumv=sumv+vpsipred(j,ndih_constr)
406 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
408 phibound(1,ndih_constr)=phihel*deg2rad
409 phibound(2,ndih_constr)=phibet*deg2rad
410 sdihed(1,ndih_constr)=sigmahel*deg2rad
411 sdihed(2,ndih_constr)=sigmabet*deg2rad
415 & 'There are',ndih_constr,
416 & ' bimodal restraints on gamma angles.'
418 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
419 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
420 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
421 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
422 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
423 & (vpsipred(j,i),j=1,3)
426 endif ! endif ndif_constr.gt.0
427 endif ! with_dihed_constr
428 if (with_theta_constr) then
429 C with_theta_constr is keyword allowing for occurance of theta constrains
430 read (inp,*) ntheta_constr
431 C ntheta_constr is the number of theta constrains
432 if (ntheta_constr.gt.0) then
434 read (inp,*) (itheta_constr(i),theta_constr0(i),
435 & theta_drange(i),for_thet_constr(i),
437 C the above code reads from 1 to ntheta_constr
438 C itheta_constr(i) residue i for which is theta_constr
439 C theta_constr0 the global minimum value
440 C theta_drange is range for which there is no energy penalty
441 C for_thet_constr is the force constant for quartic energy penalty
443 C if(me.eq.king.or..not.out1file)then
445 & 'There are',ntheta_constr,' constraints on phi angles.'
447 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
453 theta_constr0(i)=deg2rad*theta_constr0(i)
454 theta_drange(i)=deg2rad*theta_drange(i)
456 C if(me.eq.king.or..not.out1file)
457 C & write (iout,*) 'FTORS',ftors
458 C do i=1,ntheta_constr
459 C ii = itheta_constr(i)
460 C thetabound(1,ii) = phi0(i)-drange(i)
461 C thetabound(2,ii) = phi0(i)+drange(i)
463 endif ! ntheta_constr.gt.0
464 endif! with_theta_constr
467 c read(inp,'(a)') pdbfile
468 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
469 c open(ipdbin,file=pdbfile,status='old',err=33)
471 c 33 write (iout,'(a)') 'Error opening PDB file.'
474 c print *,'Begin reading pdb data'
476 c print *,'Finished reading pdb data'
477 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
479 c itype_pdb(i)=itype(i)
482 c write (iout,'(a,i3)') 'nsup=',nsup
484 c if (nsup.le.(nct-nnt+1)) then
485 c do i=0,nct-nnt+1-nsup
486 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
492 c & 'Error - sequences to be superposed do not match.'
495 c do i=0,nsup-(nct-nnt+1)
496 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
498 c nstart_sup=nstart_sup+i
504 c & 'Error - sequences to be superposed do not match.'
507 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
508 c & ' nstart_seq=',nstart_seq
512 write (iout,*) "molread: REFSTR",refstr
514 if (.not.pdbref) then
515 call read_angles(inp,*38)
517 38 write (iout,'(a)') 'Error reading reference structure.'
519 call mp_stopall(Error_Msg)
521 stop 'Error reading reference structure'
534 c call contact(.true.,ncont_ref,icont_ref)
537 C write (iout,'(/a,i3,a)')
538 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
539 write (iout,'(20i4)') (iss(i),i=1,ns)
541 write(iout,*)"Running with dynamic disulfide-bond formation"
543 write (iout,'(/a/)') 'Pre-formed links are:'
549 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
550 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
556 if (ns.gt.0.and.dyn_ss) then
560 forcon(i-nss)=forcon(i)
567 dyn_ss_mask(iss(i))=.true.
570 c Read distance restraints
571 if (constr_dist.gt.0) then
572 call read_dist_constr
577 c-----------------------------------------------------------------------------
578 logical function seq_comp(itypea,itypeb,length)
580 integer length,itypea(length),itypeb(length)
583 if (itypea(i).ne.itypeb(i)) then
591 c-----------------------------------------------------------------------------
592 subroutine read_bridge
593 C Read information about disulfide bridges.
596 include 'COMMON.IOUNITS'
599 include 'COMMON.INTERACT'
600 include 'COMMON.LOCAL'
601 include 'COMMON.NAMES'
602 include 'COMMON.CHAIN'
603 include 'COMMON.FFIELD'
604 include 'COMMON.SBRIDGE'
605 include 'COMMON.HEADER'
606 include 'COMMON.CONTROL'
607 include 'COMMON.TIME1'
609 include 'COMMON.INFO'
612 C Read bridging residues.
613 read (inp,*) ns,(iss(i),i=1,ns)
615 C Check whether the specified bridging residues are cystines.
617 if (iabs(itype(iss(i))).ne.1) then
618 write (iout,'(2a,i3,a)')
619 & 'Do you REALLY think that the residue ',
620 & restyp(itype(iss(i))),i,
621 & ' can form a disulfide bridge?!!!'
622 write (*,'(2a,i3,a)')
623 & 'Do you REALLY think that the residue ',
624 & restyp(itype(iss(i))),i,
625 & ' can form a disulfide bridge?!!!'
627 call mp_stopall(error_msg)
633 C Read preformed bridges.
635 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
638 C Check if the residues involved in bridges are in the specified list of
642 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
643 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
644 write (iout,'(a,i3,a)') 'Disulfide pair',i,
645 & ' contains residues present in other pairs.'
646 write (*,'(a,i3,a)') 'Disulfide pair',i,
647 & ' contains residues present in other pairs.'
649 call mp_stopall(error_msg)
656 if (ihpb(i).eq.iss(j)) goto 10
658 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
661 if (jhpb(i).eq.iss(j)) goto 20
663 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
676 c----------------------------------------------------------------------------
677 subroutine read_angles(kanal,*)
682 include 'COMMON.CHAIN'
683 include 'COMMON.IOUNITS'
685 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
686 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
687 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
688 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
690 theta(i)=deg2rad*theta(i)
691 phi(i)=deg2rad*phi(i)
692 alph(i)=deg2rad*alph(i)
693 omeg(i)=deg2rad*omeg(i)
698 c----------------------------------------------------------------------------
699 subroutine reada(rekord,lancuch,wartosc,default)
701 character*(*) rekord,lancuch
702 double precision wartosc,default
705 iread=index(rekord,lancuch)
710 iread=iread+ilen(lancuch)+1
711 read (rekord(iread:),*) wartosc
714 c----------------------------------------------------------------------------
715 subroutine multreada(rekord,lancuch,tablica,dim,default)
718 double precision tablica(dim),default
719 character*(*) rekord,lancuch
725 iread=index(rekord,lancuch)
726 if (iread.eq.0) return
727 iread=iread+ilen(lancuch)+1
728 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
731 c----------------------------------------------------------------------------
732 subroutine readi(rekord,lancuch,wartosc,default)
734 character*(*) rekord,lancuch
735 integer wartosc,default
738 iread=index(rekord,lancuch)
743 iread=iread+ilen(lancuch)+1
744 read (rekord(iread:),*) wartosc
747 C----------------------------------------------------------------------
748 subroutine multreadi(rekord,lancuch,tablica,dim,default)
751 integer tablica(dim),default
752 character*(*) rekord,lancuch
759 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
760 if (iread.eq.0) return
761 iread=iread+ilen(lancuch)+1
762 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
766 c----------------------------------------------------------------------------
767 subroutine card_concat(card)
769 include 'COMMON.IOUNITS'
771 character*80 karta,ucase
773 read (inp,'(a)') karta
776 do while (karta(80:80).eq.'&')
777 card=card(:ilen(card)+1)//karta(:79)
778 read (inp,'(a)') karta
781 card=card(:ilen(card)+1)//karta
784 c----------------------------------------------------------------------------
793 include 'COMMON.IOUNITS'
794 include 'COMMON.CONTROL'
795 integer lenpre,lenpot,ilen
797 character*16 cformat,cprint
799 integer lenint,lenout
800 call getenv('INPUT',prefix)
801 call getenv('OUTPUT',prefout)
802 call getenv('INTIN',prefintin)
803 call getenv('COORD',cformat)
804 call getenv('PRINTCOOR',cprint)
805 call getenv('SCRATCHDIR',scratchdir)
808 if (index(ucase(cformat),'CX').gt.0) then
815 lenint=ilen(prefintin)
816 C Get the names and open the input files
817 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
819 write (liczba,'(bz,i3.3)') me
820 outname=prefout(:lenout)//'_clust.out_'//liczba
822 outname=prefout(:lenout)//'_clust.out'
825 intinname=prefintin(:lenint)//'.bx'
826 else if (from_cx) then
827 intinname=prefintin(:lenint)//'.cx'
829 intinname=prefintin(:lenint)//'.int'
831 rmsname=prefintin(:lenint)//'.rms'
832 open (jplot,file=prefout(:ilen(prefout))//'.tex',
834 open (jrms,file=rmsname,status='unknown')
835 open(iout,file=outname,status='unknown')
836 C Get parameter filenames and open the parameter files.
837 call getenv('BONDPAR',bondname)
838 open (ibond,file=bondname,status='old')
839 call getenv('THETPAR',thetname)
840 open (ithep,file=thetname,status='old')
841 call getenv('ROTPAR',rotname)
842 open (irotam,file=rotname,status='old')
843 call getenv('TORPAR',torname)
844 open (itorp,file=torname,status='old')
846 call getenv('TORDPAR',tordname)
847 open (itordp,file=tordname,status='old')
849 call getenv('FOURIER',fouriername)
850 open (ifourier,file=fouriername,status='old')
851 call getenv('ELEPAR',elename)
852 open (ielep,file=elename,status='old')
853 call getenv('SIDEPAR',sidename)
854 open (isidep,file=sidename,status='old')
855 call getenv('SIDEP',sidepname)
856 open (isidep1,file=sidepname,status="old")
857 call getenv('SCCORPAR',sccorname)
858 open (isccor,file=sccorname,status="old")
859 call getenv('LIPTRANPAR',liptranname)
860 open (iliptranpar,file=liptranname,status='old')
863 C 8/9/01 In the newest version SCp interaction constants are read from a file
864 C Use -DOLDSCP to use hard-coded constants instead.
866 call getenv('SCPPAR',scpname)
867 open (iscpp,file=scpname,status='old')
871 c--------------------------------------------------------------------------
872 subroutine read_dist_constr
873 implicit real*8 (a-h,o-z)
878 include 'COMMON.CONTROL'
879 include 'COMMON.CHAIN'
880 include 'COMMON.IOUNITS'
881 include 'COMMON.SBRIDGE'
882 integer ifrag_(2,100),ipair_(2,100)
883 double precision wfrag_(100),wpair_(100)
884 character*500 controlcard
886 logical lprn /.true./
887 write (iout,*) "Calling read_dist_constr"
888 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
890 write(iout,*) "TU sie wywalam?"
891 call card_concat(controlcard)
892 write (iout,*) controlcard
894 call readi(controlcard,"NFRAG",nfrag_,0)
895 call readi(controlcard,"NPAIR",npair_,0)
896 call readi(controlcard,"NDIST",ndist_,0)
897 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
898 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
899 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
900 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
901 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
902 normalize = index(controlcard,"NORMALIZE").gt.0
903 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
904 write (iout,*) "IFRAG"
906 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
908 write (iout,*) "IPAIR"
910 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
914 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
915 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
916 & ifrag_(2,i)=nstart_sup+nsup-1
917 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
919 if (wfrag_(i).gt.0.0d0) then
920 do j=ifrag_(1,i),ifrag_(2,i)-1
922 write (iout,*) "j",j," k",k
924 if (constr_dist.eq.1) then
929 forcon(nhpb)=wfrag_(i)
930 else if (constr_dist.eq.2) then
931 if (ddjk.le.dist_cut) then
936 forcon(nhpb)=wfrag_(i)
943 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
946 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
947 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
953 if (wpair_(i).gt.0.0d0) then
961 do j=ifrag_(1,ii),ifrag_(2,ii)
962 do k=ifrag_(1,jj),ifrag_(2,jj)
966 forcon(nhpb)=wpair_(i)
968 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
969 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
974 write (iout,*) "Distance restraints as read from input"
976 if (constr_dist.eq.11) then
977 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
978 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
979 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
980 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
982 write (iout,'(a,4i5,2f8.2,2f10.5)') "+dist.restr ",
983 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
984 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb)
985 if (ibecarb(i).gt.0) then
991 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
992 & ibecarb(i),forcon(nhpb+1)
993 if (forcon(nhpb+1).gt.0.0d0) then
995 if (ibecarb(i).gt.0) then
999 if (dhpb(nhpb).eq.0.0d0)
1000 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1002 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1003 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1005 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1006 C if (forcon(nhpb+1).gt.0.0d0) then
1008 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1010 if (constr_dist.eq.11 .and. normalize) then
1011 fordepthmax=fordepth(1)
1013 if (fordepth(i).gt.fordepthmax) fordepthmax=fordepth(i)
1016 fordepth(i)=fordepth(i)/fordepthmax
1018 write (iout,'(/a/4a5,2a8,2a10)')
1019 & "Normalized Lorenzian-like distance restraints",
1020 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V"
1022 write (iout,'(4i5,2f8.2,2f10.5)')i,ihpb(i),jhpb(i),ibecarb(i),
1023 & dhpb(i),dhpb1(i),forcon(i),fordepth(i)