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
63 pdbref=(index(controlcard,'PDBREF').gt.0)
64 if (index(controlcard,"CASC").gt.0) then
66 else if (index(controlcard,"SCONLY").gt.0) then
72 C VSolvSphere the volume of solving sphere
74 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
75 C there will be no distinction between proline peptide group and normal peptide
76 C group in case of shielding parameters
77 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
78 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
79 write (iout,*) VSolvSphere,VSolvSphere_div
80 C long axis of side chain
82 long_r_sidechain(i)=vbldsc0(1,i)
83 short_r_sidechain(i)=sigma0(i)
87 call readi(controlcard,'PDBOUT',outpdb,0)
88 call readi(controlcard,'MOL2OUT',outmol2,0)
89 refstr=(index(controlcard,'REFSTR').gt.0)
90 pdbref=(index(controlcard,'PDBREF').gt.0)
91 refstr = refstr .or. pdbref
92 write (iout,*) "REFSTR",refstr," PDBREF",pdbref
93 iscode=index(controlcard,'ONE_LETTER')
94 tree=(index(controlcard,'MAKE_TREE').gt.0)
95 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
96 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
97 write (iout,*) "with_dihed_constr ",with_dihed_constr,
98 & " CONSTR_DIST",constr_dist
99 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
100 write (iout,*) "with_theta_constr ",with_theta_constr
102 min_var=(index(controlcard,'MINVAR').gt.0)
103 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
104 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
105 print_fittest=(index(controlcard,'PRINT_FITTEST').gt.0)
106 call readi(controlcard,'NCUT',ncut,0)
108 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
111 call readi(controlcard,'NCLUST',nclust,5)
113 call readi(controlcard,'SYM',symetr,1)
114 write (iout,*) 'sym', symetr
115 call readi(controlcard,'NSTART',nstart,0)
116 call readi(controlcard,'NEND',nend,0)
117 call reada(controlcard,'ECUT',ecut,10.0d0)
118 call reada(controlcard,'PROB',prob_limit,0.99d0)
119 write (iout,*) "Probability limit",prob_limit
120 lgrp=(index(controlcard,'LGRP').gt.0)
121 caonly=(index(controlcard,'CA_ONLY').gt.0)
122 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
123 call readi(controlcard,'IOPT',iopt,2)
124 lefree = index(controlcard,"EFREE").gt.0
125 call readi(controlcard,'NTEMP',nT,1)
126 write (iout,*) "nT",nT
127 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
128 write (iout,*) "nT",nT
129 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
131 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
133 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
134 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
135 lprint_int=index(controlcard,"PRINT_INT") .gt.0
136 call readi(controlcard,'NSAXS',nsaxs,0)
137 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
138 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
139 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
140 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
141 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
145 c--------------------------------------------------------------------------
148 C Read molecular data.
152 include 'COMMON.IOUNITS'
155 include 'COMMON.INTERACT'
156 include 'COMMON.LOCAL'
157 include 'COMMON.NAMES'
158 include 'COMMON.CHAIN'
159 include 'COMMON.FFIELD'
160 include 'COMMON.SBRIDGE'
161 include 'COMMON.HEADER'
162 include 'COMMON.CONTROL'
163 include 'COMMON.CONTACTS'
164 include 'COMMON.TIME1'
165 include 'COMMON.TORCNSTR'
166 include 'COMMON.SHIELD'
168 include 'COMMON.INFO'
170 character*4 sequence(maxres)
171 character*800 weightcard,controlcard
173 double precision x(maxvar)
174 double precision phihel,phibet,sigmahel,sigmabet,sumv,
176 integer itype_pdb(maxres)
178 integer i,j,kkk,i1,i2,it1,it2,tperm,ii,iperm
182 C Read weights of the subsequent energy terms.
183 call card_concat(weightcard)
184 call reada(weightcard,'WSC',wsc,1.0d0)
185 call reada(weightcard,'WLONG',wsc,wsc)
186 call reada(weightcard,'WSCP',wscp,1.0d0)
187 call reada(weightcard,'WELEC',welec,1.0D0)
188 call reada(weightcard,'WVDWPP',wvdwpp,welec)
189 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
190 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
191 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
192 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
193 call reada(weightcard,'WTURN3',wturn3,1.0D0)
194 call reada(weightcard,'WTURN4',wturn4,1.0D0)
195 call reada(weightcard,'WTURN6',wturn6,1.0D0)
196 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
197 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
198 call reada(weightcard,'WBOND',wbond,1.0D0)
199 call reada(weightcard,'WTOR',wtor,1.0D0)
200 call reada(weightcard,'WTORD',wtor_d,1.0D0)
201 call reada(weightcard,'WANG',wang,1.0D0)
202 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
203 call reada(weightcard,'WSAXS',wsaxs,0.0D0)
204 call reada(weightcard,'SCAL14',scal14,0.4D0)
205 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
206 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
207 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
208 if (index(weightcard,'SOFT').gt.0) ipot=6
209 call reada(weightcard,"D0CM",d0cm,3.78d0)
210 call reada(weightcard,"AKCM",akcm,15.1d0)
211 call reada(weightcard,"AKTH",akth,11.0d0)
212 call reada(weightcard,"AKCT",akct,12.0d0)
213 call reada(weightcard,"V1SS",v1ss,-1.08d0)
214 call reada(weightcard,"V2SS",v2ss,7.61d0)
215 call reada(weightcard,"V3SS",v3ss,13.7d0)
216 call reada(weightcard,"EBR",ebr,-5.50D0)
217 call reada(weightcard,'WSHIELD',wshield,1.0d0)
218 write(iout,*) 'WSHIELD',wshield
219 call reada(weightcard,'WLT',wliptran,0.0D0)
220 call reada(weightcard,"ATRISS",atriss,0.301D0)
221 call reada(weightcard,"BTRISS",btriss,0.021D0)
222 call reada(weightcard,"CTRISS",ctriss,1.001D0)
223 call reada(weightcard,"DTRISS",dtriss,1.001D0)
224 write (iout,*) "ATRISS=", atriss
225 write (iout,*) "BTRISS=", btriss
226 write (iout,*) "CTRISS=", ctriss
227 write (iout,*) "DTRISS=", dtriss
228 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
230 dyn_ss_mask(i)=.false.
234 dyn_ssbond_ij(i,j)=1.0d300
237 call reada(weightcard,"HT",Ht,0.0D0)
239 ss_depth=ebr/wsc-0.25*eps(1,1)
240 Ht=Ht/wsc-0.25*eps(1,1)
241 akcm=akcm*wstrain/wsc
242 akth=akth*wstrain/wsc
243 akct=akct*wstrain/wsc
244 v1ss=v1ss*wstrain/wsc
245 v2ss=v2ss*wstrain/wsc
246 v3ss=v3ss*wstrain/wsc
248 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
250 write (iout,'(/a)') "Disulfide bridge parameters:"
251 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
252 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
253 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
254 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
257 C 12/1/95 Added weight for the multi-body term WCORR
258 call reada(weightcard,'WCORRH',wcorr,1.0D0)
259 if (wcorr4.gt.0.0d0) wcorr=wcorr4
279 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
280 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
281 & wturn4,wturn6,wsccor
282 10 format (/'Energy-term weights (unscaled):'//
283 & 'WSCC= ',f10.6,' (SC-SC)'/
284 & 'WSCP= ',f10.6,' (SC-p)'/
285 & 'WELEC= ',f10.6,' (p-p electr)'/
286 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
287 & 'WBOND= ',f10.6,' (stretching)'/
288 & 'WANG= ',f10.6,' (bending)'/
289 & 'WSCLOC= ',f10.6,' (SC local)'/
290 & 'WTOR= ',f10.6,' (torsional)'/
291 & 'WTORD= ',f10.6,' (double torsional)'/
292 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
293 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
294 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
295 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
296 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
297 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
298 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
299 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
300 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
302 if (wcorr4.gt.0.0d0) then
303 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
304 & 'between contact pairs of peptide groups'
305 write (iout,'(2(a,f5.3/))')
306 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
307 & 'Range of quenching the correlation terms:',2*delt_corr
308 else if (wcorr.gt.0.0d0) then
309 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
310 & 'between contact pairs of peptide groups'
312 write (iout,'(a,f8.3)')
313 & 'Scaling factor of 1,4 SC-p interactions:',scal14
314 write (iout,'(a,f8.3)')
315 & 'General scaling factor of SC-p interactions:',scalscp
316 r0_corr=cutoff_corr-delt_corr
318 aad(i,1)=scalscp*aad(i,1)
319 aad(i,2)=scalscp*aad(i,2)
320 bad(i,1)=scalscp*bad(i,1)
321 bad(i,2)=scalscp*bad(i,2)
325 c print *,'indpdb=',indpdb,' pdbref=',pdbref
327 C Read sequence if not taken from the pdb file.
328 if (iscode.gt.0) then
329 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
331 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
333 C Convert sequence to numeric code
335 itype(i)=rescode(i,sequence(i),iscode)
338 c print '(20i4)',(itype(i),i=1,nres)
342 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
344 if (itype(i).eq.ntyp1) then
348 else if (iabs(itype(i+1)).ne.20) then
350 else if (iabs(itype(i)).ne.20) then
357 write (iout,*) "ITEL"
359 write (iout,*) i,itype(i),itel(i)
362 c print *,'Call Read_Bridge.'
364 C this fragment reads diheadral constrains
367 c print *,'NNT=',NNT,' NCT=',NCT
368 call seq2chains(nres,itype,nchain,chain_length,chain_border,
370 write(iout,*) "nres",nres," nchain",nchain
372 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
375 call chain_symmetry(nchain,nres,itype,chain_border,
376 & chain_length,npermchain,tabpermchain)
378 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
381 write(iout,*) "residue permutations"
383 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
385 if (itype(1).eq.ntyp1) nnt=2
386 if (itype(nres).eq.ntyp1) nct=nct-1
387 if (nstart.lt.nnt) nstart=nnt
388 if (nend.gt.nct .or. nend.eq.0) nend=nct
389 write (iout,*) "nstart",nstart," nend",nend
391 if (with_dihed_constr) then
393 read (inp,*) ndih_constr
394 if (ndih_constr.gt.0) then
397 C write (iout,*) 'FTORS',ftors
398 C ftors is the force constant for torsional quartic constrains
399 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
402 & 'There are',ndih_constr,' constraints on phi angles.'
404 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
408 phi0(i)=deg2rad*phi0(i)
409 drange(i)=deg2rad*drange(i)
411 else if (ndih_constr.lt.0) then
413 call card_concat(controlcard)
414 call reada(controlcard,"PHIHEL",phihel,50.0D0)
415 call reada(controlcard,"PHIBET",phibet,180.0D0)
416 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
417 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
418 call reada(controlcard,"WDIHC",wdihc,0.591d0)
419 write (iout,*) "Weight of the dihedral restraint term",wdihc
420 read(inp,'(9x,3f7.3)')
421 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
422 write (iout,*) "The secprob array"
424 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
428 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
429 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
430 ndih_constr=ndih_constr+1
431 idih_constr(ndih_constr)=i
434 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
435 sumv=sumv+vpsipred(j,ndih_constr)
438 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
440 phibound(1,ndih_constr)=phihel*deg2rad
441 phibound(2,ndih_constr)=phibet*deg2rad
442 sdihed(1,ndih_constr)=sigmahel*deg2rad
443 sdihed(2,ndih_constr)=sigmabet*deg2rad
447 & 'There are',ndih_constr,
448 & ' bimodal restraints on gamma angles.'
450 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
451 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
452 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
453 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
454 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
455 & (vpsipred(j,i),j=1,3)
458 endif ! endif ndif_constr.gt.0
459 endif ! with_dihed_constr
460 if (with_theta_constr) then
461 C with_theta_constr is keyword allowing for occurance of theta constrains
462 read (inp,*) ntheta_constr
463 C ntheta_constr is the number of theta constrains
464 if (ntheta_constr.gt.0) then
466 read (inp,*) (itheta_constr(i),theta_constr0(i),
467 & theta_drange(i),for_thet_constr(i),
469 C the above code reads from 1 to ntheta_constr
470 C itheta_constr(i) residue i for which is theta_constr
471 C theta_constr0 the global minimum value
472 C theta_drange is range for which there is no energy penalty
473 C for_thet_constr is the force constant for quartic energy penalty
476 & 'There are',ntheta_constr,' constraints on phi angles.'
478 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
484 theta_constr0(i)=deg2rad*theta_constr0(i)
485 theta_drange(i)=deg2rad*theta_drange(i)
487 C if(me.eq.king.or..not.out1file)
488 C & write (iout,*) 'FTORS',ftors
489 C do i=1,ntheta_constr
490 C ii = itheta_constr(i)
491 C thetabound(1,ii) = phi0(i)-drange(i)
492 C thetabound(2,ii) = phi0(i)+drange(i)
494 endif ! ntheta_constr.gt.0
495 endif! with_theta_constr
496 write (iout,*) "nstart",nstart," nend",nend
497 write (iout,*) "calling read_saxs_consrtr",nsaxs
498 if (nsaxs.gt.0) call read_saxs_constr
501 c read(inp,'(a)') pdbfile
502 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
503 c open(ipdbin,file=pdbfile,status='old',err=33)
505 c 33 write (iout,'(a)') 'Error opening PDB file.'
508 c print *,'Begin reading pdb data'
510 c print *,'Finished reading pdb data'
511 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
513 c itype_pdb(i)=itype(i)
516 c write (iout,'(a,i3)') 'nsup=',nsup
518 c if (nsup.le.(nct-nnt+1)) then
519 c do i=0,nct-nnt+1-nsup
520 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
526 c & 'Error - sequences to be superposed do not match.'
529 c do i=0,nsup-(nct-nnt+1)
530 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
532 c nstart_sup=nstart_sup+i
538 c & 'Error - sequences to be superposed do not match.'
541 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
542 c & ' nstart_seq=',nstart_seq
546 write (iout,*) "molread: REFSTR",refstr
548 if (.not.pdbref) then
549 call read_angles(inp,*38)
551 38 write (iout,'(a)') 'Error reading reference structure.'
553 call mp_stopall(Error_Msg)
555 stop 'Error reading reference structure'
567 c call contact(.true.,ncont_ref,icont_ref)
570 C write (iout,'(/a,i3,a)')
571 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
572 write (iout,'(20i4)') (iss(i),i=1,ns)
574 write(iout,*)"Running with dynamic disulfide-bond formation"
576 write (iout,'(/a/)') 'Pre-formed links are:'
582 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
583 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
589 if (ns.gt.0.and.dyn_ss) then
593 forcon(i-nss)=forcon(i)
600 dyn_ss_mask(iss(i))=.true.
603 c Read distance restraints
604 if (constr_dist.gt.0) then
605 call read_dist_constr
610 c-----------------------------------------------------------------------------
611 logical function seq_comp(itypea,itypeb,length)
613 integer length,itypea(length),itypeb(length)
616 if (itypea(i).ne.itypeb(i)) then
624 c-----------------------------------------------------------------------------
625 subroutine read_bridge
626 C Read information about disulfide bridges.
629 include 'COMMON.IOUNITS'
632 include 'COMMON.INTERACT'
633 include 'COMMON.LOCAL'
634 include 'COMMON.NAMES'
635 include 'COMMON.CHAIN'
636 include 'COMMON.FFIELD'
637 include 'COMMON.SBRIDGE'
638 include 'COMMON.HEADER'
639 include 'COMMON.CONTROL'
640 include 'COMMON.TIME1'
642 include 'COMMON.INFO'
645 C Read bridging residues.
646 read (inp,*) ns,(iss(i),i=1,ns)
648 C Check whether the specified bridging residues are cystines.
650 if (itype(iss(i)).ne.1) then
651 write (iout,'(2a,i3,a)')
652 & 'Do you REALLY think that the residue ',
653 & restyp(itype(iss(i))),i,
654 & ' can form a disulfide bridge?!!!'
655 write (*,'(2a,i3,a)')
656 & 'Do you REALLY think that the residue ',
657 & restyp(itype(iss(i))),i,
658 & ' can form a disulfide bridge?!!!'
660 call mp_stopall(error_msg)
666 C Read preformed bridges.
668 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
671 C Check if the residues involved in bridges are in the specified list of
675 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
676 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
677 write (iout,'(a,i3,a)') 'Disulfide pair',i,
678 & ' contains residues present in other pairs.'
679 write (*,'(a,i3,a)') 'Disulfide pair',i,
680 & ' contains residues present in other pairs.'
682 call mp_stopall(error_msg)
689 if (ihpb(i).eq.iss(j)) goto 10
691 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
694 if (jhpb(i).eq.iss(j)) goto 20
696 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
709 c----------------------------------------------------------------------------
710 subroutine read_angles(kanal,*)
715 include 'COMMON.CHAIN'
716 include 'COMMON.IOUNITS'
718 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
719 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
720 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
721 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
723 theta(i)=deg2rad*theta(i)
724 phi(i)=deg2rad*phi(i)
725 alph(i)=deg2rad*alph(i)
726 omeg(i)=deg2rad*omeg(i)
731 c----------------------------------------------------------------------------
732 subroutine reada(rekord,lancuch,wartosc,default)
734 character*(*) rekord,lancuch
735 double precision wartosc,default
738 iread=index(rekord,lancuch)
743 iread=iread+ilen(lancuch)+1
744 read (rekord(iread:),*) wartosc
747 c----------------------------------------------------------------------------
748 subroutine multreada(rekord,lancuch,tablica,dim,default)
751 double precision tablica(dim),default
752 character*(*) rekord,lancuch
758 iread=index(rekord,lancuch)
759 if (iread.eq.0) return
760 iread=iread+ilen(lancuch)+1
761 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
764 c----------------------------------------------------------------------------
765 subroutine readi(rekord,lancuch,wartosc,default)
767 character*(*) rekord,lancuch
768 integer wartosc,default
771 iread=index(rekord,lancuch)
776 iread=iread+ilen(lancuch)+1
777 read (rekord(iread:),*) wartosc
780 C----------------------------------------------------------------------
781 subroutine multreadi(rekord,lancuch,tablica,dim,default)
784 integer tablica(dim),default
785 character*(*) rekord,lancuch
792 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
793 if (iread.eq.0) return
794 iread=iread+ilen(lancuch)+1
795 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
799 c----------------------------------------------------------------------------
800 subroutine card_concat(card)
802 include 'COMMON.IOUNITS'
804 character*80 karta,ucase
806 read (inp,'(a)') karta
809 do while (karta(80:80).eq.'&')
810 card=card(:ilen(card)+1)//karta(:79)
811 read (inp,'(a)') karta
814 card=card(:ilen(card)+1)//karta
817 c----------------------------------------------------------------------------
826 include 'COMMON.IOUNITS'
827 include 'COMMON.CONTROL'
828 integer lenpre,lenpot,ilen
830 character*16 cformat,cprint
832 integer lenint,lenout
833 call getenv('INPUT',prefix)
834 call getenv('OUTPUT',prefout)
835 call getenv('INTIN',prefintin)
836 call getenv('COORD',cformat)
837 call getenv('PRINTCOOR',cprint)
838 call getenv('SCRATCHDIR',scratchdir)
841 if (index(ucase(cformat),'CX').gt.0) then
848 lenint=ilen(prefintin)
849 C Get the names and open the input files
850 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
852 write (liczba,'(bz,i3.3)') me
853 outname=prefout(:lenout)//'_clust.out_'//liczba
855 outname=prefout(:lenout)//'_clust.out'
858 intinname=prefintin(:lenint)//'.bx'
859 else if (from_cx) then
860 intinname=prefintin(:lenint)//'.cx'
862 intinname=prefintin(:lenint)//'.int'
864 rmsname=prefintin(:lenint)//'.rms'
865 open (jplot,file=prefout(:ilen(prefout))//'.tex',
867 open (jrms,file=rmsname,status='unknown')
868 open(iout,file=outname,status='unknown')
869 C Get parameter filenames and open the parameter files.
870 call getenv('BONDPAR',bondname)
871 open (ibond,file=bondname,status='old')
872 call getenv('THETPAR',thetname)
873 open (ithep,file=thetname,status='old')
874 call getenv('ROTPAR',rotname)
875 open (irotam,file=rotname,status='old')
876 call getenv('TORPAR',torname)
877 open (itorp,file=torname,status='old')
879 call getenv('TORDPAR',tordname)
880 open (itordp,file=tordname,status='old')
882 call getenv('FOURIER',fouriername)
883 open (ifourier,file=fouriername,status='old')
884 call getenv('ELEPAR',elename)
885 open (ielep,file=elename,status='old')
886 call getenv('SIDEPAR',sidename)
887 open (isidep,file=sidename,status='old')
888 call getenv('SIDEP',sidepname)
889 open (isidep1,file=sidepname,status="old")
890 call getenv('SCCORPAR',sccorname)
891 open (isccor,file=sccorname,status="old")
892 call getenv('LIPTRANPAR',liptranname)
893 open (iliptranpar,file=liptranname,status='old')
896 C 8/9/01 In the newest version SCp interaction constants are read from a file
897 C Use -DOLDSCP to use hard-coded constants instead.
899 call getenv('SCPPAR',scpname)
900 open (iscpp,file=scpname,status='old')
904 c--------------------------------------------------------------------------
905 subroutine read_dist_constr
906 implicit real*8 (a-h,o-z)
911 include 'COMMON.CONTROL'
912 include 'COMMON.CHAIN'
913 include 'COMMON.IOUNITS'
914 include 'COMMON.SBRIDGE'
915 include 'COMMON.INTERACT'
916 integer ifrag_(2,100),ipair_(2,100)
917 double precision wfrag_(100),wpair_(100)
918 character*500 controlcard
919 logical lprn /.true./
920 logical normalize,next
922 double precision scal_bfac
923 double precision xlink(4,0:4) /
925 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
926 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
927 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
928 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
929 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
930 write (iout,*) "Calling read_dist_constr"
931 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
937 call card_concat(controlcard)
938 next = index(controlcard,"NEXT").gt.0
939 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
940 write (iout,*) "restr_type",restr_type
941 call readi(controlcard,"NFRAG",nfrag_,0)
942 call readi(controlcard,"NFRAG",nfrag_,0)
943 call readi(controlcard,"NPAIR",npair_,0)
944 call readi(controlcard,"NDIST",ndist_,0)
945 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
946 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
947 if (restr_type.eq.10)
948 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
949 if (restr_type.eq.12)
950 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
951 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
952 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
953 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
954 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
955 normalize = index(controlcard,"NORMALIZE").gt.0
956 write (iout,*) "WBOLTZD",wboltzd
957 write (iout,*) "SCAL_PEAK",scal_peak
958 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
959 write (iout,*) "IFRAG"
961 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
963 write (iout,*) "IPAIR"
965 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
967 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
969 read(inp,'(a)') pdbfile
971 & "Distance restraints will be constructed from structure ",pdbfile
972 open(ipdbin,file=pdbfile,status='old',err=11)
978 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
979 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
980 & ifrag_(2,i)=nstart_sup+nsup-1
981 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
983 if (wfrag_(i).eq.0.0d0) cycle
984 do j=ifrag_(1,i),ifrag_(2,i)-1
986 c write (iout,*) "j",j," k",k
988 if (restr_type.eq.1) then
994 forcon(nhpb)=wfrag_(i)
995 else if (constr_dist.eq.2) then
996 if (ddjk.le.dist_cut) then
1002 forcon(nhpb)=wfrag_(i)
1004 else if (restr_type.eq.3) then
1010 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1012 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
1013 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1018 if (wpair_(i).eq.0.0d0) cycle
1026 do j=ifrag_(1,ii),ifrag_(2,ii)
1027 do k=ifrag_(1,jj),ifrag_(2,jj)
1029 if (restr_type.eq.1) then
1035 forcon(nhpb)=wpair_(i)
1036 else if (constr_dist.eq.2) then
1037 if (ddjk.le.dist_cut) then
1043 forcon(nhpb)=wpair_(i)
1045 else if (restr_type.eq.3) then
1051 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1053 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1054 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1060 write (iout,*) "Distance restraints as read from input"
1062 if (restr_type.eq.12) then
1063 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1064 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1065 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1066 & fordepth_peak(nhpb_peak+1),npeak
1067 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1068 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1069 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1070 c & fordepth_peak(nhpb_peak+1),npeak
1071 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1072 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1073 nhpb_peak=nhpb_peak+1
1074 irestr_type_peak(nhpb_peak)=12
1075 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1077 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1078 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1079 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1080 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1081 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1082 if (ibecarb_peak(nhpb_peak).eq.3) then
1083 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1084 else if (ibecarb_peak(nhpb_peak).eq.2) then
1085 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1086 else if (ibecarb_peak(nhpb_peak).eq.1) then
1087 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1088 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1090 else if (restr_type.eq.11) then
1091 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1092 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1093 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1094 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1096 irestr_type(nhpb)=11
1097 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1098 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1099 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1100 c if (ibecarb(nhpb).gt.0) then
1101 c ihpb(nhpb)=ihpb(nhpb)+nres
1102 c jhpb(nhpb)=jhpb(nhpb)+nres
1104 if (ibecarb(nhpb).eq.3) then
1105 ihpb(nhpb)=ihpb(nhpb)+nres
1106 else if (ibecarb(nhpb).eq.2) then
1107 ihpb(nhpb)=ihpb(nhpb)+nres
1108 else if (ibecarb(nhpb).eq.1) then
1109 ihpb(nhpb)=ihpb(nhpb)+nres
1110 jhpb(nhpb)=jhpb(nhpb)+nres
1112 else if (restr_type.eq.10) then
1113 c Cross-lonk Markov-like potential
1114 call card_concat(controlcard)
1115 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1116 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1118 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1119 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1120 if (index(controlcard,"ZL").gt.0) then
1122 else if (index(controlcard,"ADH").gt.0) then
1124 else if (index(controlcard,"PDH").gt.0) then
1126 else if (index(controlcard,"DSS").gt.0) then
1131 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1132 & xlink(1,link_type))
1133 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1134 & xlink(2,link_type))
1135 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1136 & xlink(3,link_type))
1137 call reada(controlcard,"SIGMA",forcon(nhpb+1),
1138 & xlink(4,link_type))
1139 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1140 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1141 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1142 if (forcon(nhpb+1).le.0.0d0 .or.
1143 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1145 irestr_type(nhpb)=10
1146 c if (ibecarb(nhpb).gt.0) then
1147 c ihpb(nhpb)=ihpb(nhpb)+nres
1148 c jhpb(nhpb)=jhpb(nhpb)+nres
1150 if (ibecarb(nhpb).eq.3) then
1151 jhpb(nhpb)=jhpb(nhpb)+nres
1152 else if (ibecarb(nhpb).eq.2) then
1153 ihpb(nhpb)=ihpb(nhpb)+nres
1154 else if (ibecarb(nhpb).eq.1) then
1155 ihpb(nhpb)=ihpb(nhpb)+nres
1156 jhpb(nhpb)=jhpb(nhpb)+nres
1158 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1159 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1160 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1164 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1165 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1166 if (forcon(nhpb+1).gt.0.0d0) then
1168 if (dhpb1(nhpb).eq.0.0d0) then
1173 c if (ibecarb(nhpb).gt.0) then
1174 c ihpb(nhpb)=ihpb(nhpb)+nres
1175 c jhpb(nhpb)=jhpb(nhpb)+nres
1177 if (ibecarb(nhpb).eq.3) then
1178 jhpb(nhpb)=jhpb(nhpb)+nres
1179 else if (ibecarb(nhpb).eq.2) then
1180 ihpb(nhpb)=ihpb(nhpb)+nres
1181 else if (ibecarb(nhpb).eq.1) then
1182 ihpb(nhpb)=ihpb(nhpb)+nres
1183 jhpb(nhpb)=jhpb(nhpb)+nres
1185 if (dhpb(nhpb).eq.0.0d0)
1186 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1188 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1189 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1191 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1192 C if (forcon(nhpb+1).gt.0.0d0) then
1194 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1197 if (restr_type.eq.4) then
1198 write (iout,*) "The BFAC array"
1200 write (iout,'(i5,f10.5)') i,bfac(i)
1203 if (itype(i).eq.ntyp1) cycle
1205 if (itype(j).eq.ntyp1) cycle
1206 if (itype(i).eq.10) then
1211 if (itype(j).eq.10) then
1221 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1224 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1225 ihpb(nhpb)=i+nres*ii
1226 jhpb(nhpb)=j+nres*jj
1227 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1228 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1229 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1230 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1239 if (restr_type.eq.5) then
1240 restr_on_coord=.true.
1242 if (itype(i).eq.ntyp1) cycle
1243 bfac(i)=(scal_bfac/bfac(i))**2
1252 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1253 & fordepthmax=fordepth(i)
1256 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1259 if (nhpb.gt.nss) then
1260 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1261 & "The following",nhpb-nss,
1262 & " distance restraints have been imposed:",
1263 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1266 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1267 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1274 11 write (iout,*)"read_dist_restr: error reading reference structure"
1277 c-------------------------------------------------------------------------------
1278 subroutine read_saxs_constr
1279 implicit real*8 (a-h,o-z)
1280 include 'DIMENSIONS'
1284 include 'COMMON.CONTROL'
1285 include 'COMMON.CHAIN'
1286 include 'COMMON.IOUNITS'
1287 include 'COMMON.SBRIDGE'
1288 double precision cm(3)
1290 write (iout,*) "Calling read_saxs nsaxs",nsaxs
1292 if (saxs_mode.eq.0) then
1293 c SAXS distance distribution
1295 read(inp,*) distsaxs(i),Psaxs(i)
1299 Cnorm = Cnorm + Psaxs(i)
1301 write (iout,*) "Cnorm",Cnorm
1303 Psaxs(i)=Psaxs(i)/Cnorm
1305 write (iout,*) "Normalized distance distribution from SAXS"
1307 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1311 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1313 write (iout,*) "Wsaxs0",Wsaxs0
1317 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1324 cm(j)=cm(j)+Csaxs(j,i)
1332 Csaxs(j,i)=Csaxs(j,i)-cm(j)
1335 write (iout,*) "SAXS sphere coordinates"
1337 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)