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'
21 character*320 controlcard,ucase
25 integer i,i1,i2,it1,it2
27 read (INP,'(a80)') titel
28 call card_concat(controlcard)
30 dyn_ss = index(controlcard,"DYN_SS").gt.0
31 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
32 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
33 call readi(controlcard,'TORMODE',tor_mode,0)
34 call readi(controlcard,'NRES',nres,0)
35 call readi(controlcard,'RESCALE',rescale_mode,2)
36 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
37 write (iout,*) "DISTCHAINMAX",distchainmax
38 C Reading the dimensions of box in x,y,z coordinates
39 call reada(controlcard,'BOXX',boxxsize,100.0d0)
40 call reada(controlcard,'BOXY',boxysize,100.0d0)
41 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
42 c Cutoff range for interactions
43 call reada(controlcard,"R_CUT",r_cut,25.0d0)
44 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
45 write (iout,*) "Cutoff on interactions",r_cut
46 write (iout,*) "lambda",rlamb
47 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
48 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
49 if (lipthick.gt.0.0d0) then
50 bordliptop=(boxzsize+lipthick)/2.0
51 bordlipbot=bordliptop-lipthick
53 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
54 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
55 buflipbot=bordlipbot+lipbufthick
56 bufliptop=bordliptop-lipbufthick
57 if ((lipbufthick*2.0d0).gt.lipthick)
58 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
60 write(iout,*) "bordliptop=",bordliptop
61 write(iout,*) "bordlipbot=",bordlipbot
62 write(iout,*) "bufliptop=",bufliptop
63 write(iout,*) "buflipbot=",buflipbot
65 call readi(controlcard,'SHIELD',shield_mode,0)
66 write (iout,*) "SHIELD MODE",shield_mode
67 if (shield_mode.gt.0) then
68 pdbref=(index(controlcard,'PDBREF').gt.0)
69 if (index(controlcard,"CASC").gt.0) then
71 else if (index(controlcard,"SCONLY").gt.0) then
77 C VSolvSphere the volume of solving sphere
79 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
80 C there will be no distinction between proline peptide group and normal peptide
81 C group in case of shielding parameters
82 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
83 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
84 write (iout,*) VSolvSphere,VSolvSphere_div
85 C long axis of side chain
87 long_r_sidechain(i)=vbldsc0(1,i)
88 short_r_sidechain(i)=sigma0(i)
92 call readi(controlcard,'PDBOUT',outpdb,0)
93 call readi(controlcard,'MOL2OUT',outmol2,0)
94 refstr=(index(controlcard,'REFSTR').gt.0)
95 pdbref=(index(controlcard,'PDBREF').gt.0)
96 refstr = refstr .or. pdbref
97 write (iout,*) "REFSTR",refstr," PDBREF",pdbref
98 iscode=index(controlcard,'ONE_LETTER')
99 tree=(index(controlcard,'MAKE_TREE').gt.0)
100 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
101 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
102 write (iout,*) "with_dihed_constr ",with_dihed_constr,
103 & " CONSTR_DIST",constr_dist
104 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
105 write (iout,*) "with_theta_constr ",with_theta_constr
107 min_var=(index(controlcard,'MINVAR').gt.0)
108 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
109 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
110 print_fittest=(index(controlcard,'PRINT_FITTEST').gt.0)
111 call readi(controlcard,'NCUT',ncut,0)
113 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
116 call readi(controlcard,'NCLUST',nclust,5)
118 call readi(controlcard,'SYM',symetr,1)
119 write (iout,*) 'sym', symetr
120 call readi(controlcard,'NSTART',nstart,0)
121 call readi(controlcard,'NEND',nend,0)
122 call reada(controlcard,'ECUT',ecut,10.0d0)
123 call reada(controlcard,'PROB',prob_limit,0.99d0)
124 write (iout,*) "Probability limit",prob_limit
125 lgrp=(index(controlcard,'LGRP').gt.0)
126 caonly=(index(controlcard,'CA_ONLY').gt.0)
127 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
128 call readi(controlcard,'IOPT',iopt,2)
129 lefree = index(controlcard,"EFREE").gt.0
130 call readi(controlcard,'NTEMP',nT,1)
131 write (iout,*) "nT",nT
132 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
133 write (iout,*) "nT",nT
134 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
136 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
138 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
139 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
140 lprint_int=index(controlcard,"PRINT_INT") .gt.0
141 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
142 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
143 c if (constr_homology) tole=dmax1(tole,1.5d0)
144 write (iout,*) "with_homology_constr ",with_dihed_constr,
145 & " CONSTR_HOMOLOGY",constr_homology
146 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
147 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
148 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
149 write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD
150 call readi(controlcard,'NSAXS',nsaxs,0)
151 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
152 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
153 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
154 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
155 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
159 c--------------------------------------------------------------------------
162 C Read molecular data.
166 include 'COMMON.IOUNITS'
169 include 'COMMON.INTERACT'
170 include 'COMMON.LOCAL'
171 include 'COMMON.NAMES'
172 include 'COMMON.CHAIN'
173 include 'COMMON.FFIELD'
174 include 'COMMON.SBRIDGE'
175 include 'COMMON.HEADER'
176 include 'COMMON.CONTROL'
177 include 'COMMON.CONTACTS'
178 include 'COMMON.TIME1'
179 include 'COMMON.TORCNSTR'
180 include 'COMMON.SHIELD'
181 include 'COMMON.SAXS'
183 include 'COMMON.INFO'
185 character*4 sequence(maxres)
186 character*800 weightcard,controlcard
188 double precision x(maxvar)
189 double precision phihel,phibet,sigmahel,sigmabet,sumv,
191 integer itype_pdb(maxres)
193 integer i,j,kkk,i1,i2,it1,it2,tperm,ii,iperm
197 C Read weights of the subsequent energy terms.
198 call card_concat(weightcard)
199 call reada(weightcard,'WSC',wsc,1.0d0)
200 call reada(weightcard,'WLONG',wsc,wsc)
201 call reada(weightcard,'WSCP',wscp,1.0d0)
202 call reada(weightcard,'WELEC',welec,1.0D0)
203 call reada(weightcard,'WVDWPP',wvdwpp,welec)
204 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
205 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
206 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
207 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
208 call reada(weightcard,'WTURN3',wturn3,1.0D0)
209 call reada(weightcard,'WTURN4',wturn4,1.0D0)
210 call reada(weightcard,'WTURN6',wturn6,1.0D0)
211 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
212 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
213 call reada(weightcard,'WBOND',wbond,1.0D0)
214 call reada(weightcard,'WTOR',wtor,1.0D0)
215 call reada(weightcard,'WTORD',wtor_d,1.0D0)
216 call reada(weightcard,'WANG',wang,1.0D0)
217 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
218 call reada(weightcard,'WSAXS',wsaxs,0.0D0)
219 call reada(weightcard,'SCAL14',scal14,0.4D0)
220 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
221 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
222 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
223 if (index(weightcard,'SOFT').gt.0) ipot=6
224 call reada(weightcard,"D0CM",d0cm,3.78d0)
225 call reada(weightcard,"AKCM",akcm,15.1d0)
226 call reada(weightcard,"AKTH",akth,11.0d0)
227 call reada(weightcard,"AKCT",akct,12.0d0)
228 call reada(weightcard,"V1SS",v1ss,-1.08d0)
229 call reada(weightcard,"V2SS",v2ss,7.61d0)
230 call reada(weightcard,"V3SS",v3ss,13.7d0)
231 call reada(weightcard,"EBR",ebr,-5.50D0)
232 call reada(weightcard,'WSHIELD',wshield,1.0d0)
233 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
234 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
235 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
236 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
237 call reada(weightcard,'WLT',wliptran,0.0D0)
238 call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
239 call reada(weightcard,"ATRISS",atriss,0.301D0)
240 call reada(weightcard,"BTRISS",btriss,0.021D0)
241 call reada(weightcard,"CTRISS",ctriss,1.001D0)
242 call reada(weightcard,"DTRISS",dtriss,1.001D0)
243 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
245 dyn_ss_mask(i)=.false.
249 dyn_ssbond_ij(i,j)=1.0d300
252 call reada(weightcard,"HT",Ht,0.0D0)
254 ss_depth=ebr/wsc-0.25*eps(1,1)
255 Ht=Ht/wsc-0.25*eps(1,1)
256 akcm=akcm*wstrain/wsc
257 akth=akth*wstrain/wsc
258 akct=akct*wstrain/wsc
259 v1ss=v1ss*wstrain/wsc
260 v2ss=v2ss*wstrain/wsc
261 v3ss=v3ss*wstrain/wsc
263 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
265 write (iout,'(/a)') "Disulfide bridge parameters:"
266 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
267 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
268 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
269 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
271 write (iout,*) "Parameters of the 'trisulfide' potential"
272 write (iout,*) "ATRISS=", atriss
273 write (iout,*) "BTRISS=", btriss
274 write (iout,*) "CTRISS=", ctriss
275 write (iout,*) "DTRISS=", dtriss
277 C 12/1/95 Added weight for the multi-body term WCORR
278 call reada(weightcard,'WCORRH',wcorr,1.0D0)
279 if (wcorr4.gt.0.0d0) wcorr=wcorr4
299 weights(28)=wdfa_dist
302 weights(31)=wdfa_beta
303 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
304 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
305 & wturn4,wturn6,wsccor
306 10 format (/'Energy-term weights (unscaled):'//
307 & 'WSCC= ',f10.6,' (SC-SC)'/
308 & 'WSCP= ',f10.6,' (SC-p)'/
309 & 'WELEC= ',f10.6,' (p-p electr)'/
310 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
311 & 'WBOND= ',f10.6,' (stretching)'/
312 & 'WANG= ',f10.6,' (bending)'/
313 & 'WSCLOC= ',f10.6,' (SC local)'/
314 & 'WTOR= ',f10.6,' (torsional)'/
315 & 'WTORD= ',f10.6,' (double torsional)'/
316 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
317 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
318 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
319 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
320 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
321 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
322 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
323 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
324 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
326 if (wcorr4.gt.0.0d0) then
327 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
328 & 'between contact pairs of peptide groups'
329 write (iout,'(2(a,f5.3/))')
330 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
331 & 'Range of quenching the correlation terms:',2*delt_corr
332 else if (wcorr.gt.0.0d0) then
333 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
334 & 'between contact pairs of peptide groups'
336 write (iout,'(a,f8.3)')
337 & 'Scaling factor of 1,4 SC-p interactions:',scal14
338 write (iout,'(a,f8.3)')
339 & 'General scaling factor of SC-p interactions:',scalscp
340 r0_corr=cutoff_corr-delt_corr
342 aad(i,1)=scalscp*aad(i,1)
343 aad(i,2)=scalscp*aad(i,2)
344 bad(i,1)=scalscp*bad(i,1)
345 bad(i,2)=scalscp*bad(i,2)
348 write (iout,'(/a/)') "DFA pseudopotential parameters:"
349 write (iout,'(a,f10.6,a)')
350 & "WDFAD= ",wdfa_dist," (distance)",
351 & "WDFAT= ",wdfa_tor," (backbone angles)",
352 & "WDFAN= ",wdfa_nei," (neighbors)",
353 & "WDFAB= ",wdfa_beta," (beta structure)"
356 c print *,'indpdb=',indpdb,' pdbref=',pdbref
358 C Read sequence if not taken from the pdb file.
359 if (iscode.gt.0) then
360 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
362 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
364 C Convert sequence to numeric code
366 itype(i)=rescode(i,sequence(i),iscode)
369 c print '(20i4)',(itype(i),i=1,nres)
373 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
375 if (itype(i).eq.ntyp1) then
379 else if (iabs(itype(i+1)).ne.20) then
381 else if (iabs(itype(i)).ne.20) then
388 write (iout,*) "ITEL"
390 write (iout,*) i,itype(i),itel(i)
393 c print *,'Call Read_Bridge.'
395 C this fragment reads diheadral constrains
398 c print *,'NNT=',NNT,' NCT=',NCT
399 call seq2chains(nres,itype,nchain,chain_length,chain_border,
401 write(iout,*) "nres",nres," nchain",nchain
403 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
406 call chain_symmetry(nchain,nres,itype,chain_border,
407 & chain_length,npermchain,tabpermchain)
409 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
412 write(iout,*) "residue permutations"
414 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
416 if (itype(1).eq.ntyp1) nnt=2
417 if (itype(nres).eq.ntyp1) nct=nct-1
418 if (nstart.lt.nnt) nstart=nnt
419 if (nend.gt.nct .or. nend.eq.0) nend=nct
420 write (iout,*) "nstart",nstart," nend",nend
423 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
424 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
426 print*, 'init_dfa_vars finished!'
428 print*, 'read_dfa_info finished!'
431 C If the reference structure is not read set the superposition
438 if (with_dihed_constr) then
440 read (inp,*) ndih_constr
441 if (ndih_constr.gt.0) then
444 C write (iout,*) 'FTORS',ftors
445 C ftors is the force constant for torsional quartic constrains
446 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
449 & 'There are',ndih_constr,' constraints on phi angles.'
451 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
455 phi0(i)=deg2rad*phi0(i)
456 drange(i)=deg2rad*drange(i)
458 else if (ndih_constr.lt.0) then
460 call card_concat(controlcard)
461 call reada(controlcard,"PHIHEL",phihel,50.0D0)
462 call reada(controlcard,"PHIBET",phibet,180.0D0)
463 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
464 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
465 call reada(controlcard,"WDIHC",wdihc,0.591d0)
466 write (iout,*) "Weight of the dihedral restraint term",wdihc
467 read(inp,'(9x,3f7.3)')
468 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
469 write (iout,*) "The secprob array"
471 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
475 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
476 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
477 ndih_constr=ndih_constr+1
478 idih_constr(ndih_constr)=i
481 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
482 sumv=sumv+vpsipred(j,ndih_constr)
485 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
487 phibound(1,ndih_constr)=phihel*deg2rad
488 phibound(2,ndih_constr)=phibet*deg2rad
489 sdihed(1,ndih_constr)=sigmahel*deg2rad
490 sdihed(2,ndih_constr)=sigmabet*deg2rad
494 & 'There are',ndih_constr,
495 & ' bimodal restraints on gamma angles.'
497 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
498 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
499 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
500 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
501 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
502 & (vpsipred(j,i),j=1,3)
505 endif ! endif ndif_constr.gt.0
506 endif ! with_dihed_constr
507 if (with_theta_constr) then
508 C with_theta_constr is keyword allowing for occurance of theta constrains
509 read (inp,*) ntheta_constr
510 C ntheta_constr is the number of theta constrains
511 if (ntheta_constr.gt.0) then
513 read (inp,*) (itheta_constr(i),theta_constr0(i),
514 & theta_drange(i),for_thet_constr(i),
516 C the above code reads from 1 to ntheta_constr
517 C itheta_constr(i) residue i for which is theta_constr
518 C theta_constr0 the global minimum value
519 C theta_drange is range for which there is no energy penalty
520 C for_thet_constr is the force constant for quartic energy penalty
523 & 'There are',ntheta_constr,' constraints on phi angles.'
525 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
531 theta_constr0(i)=deg2rad*theta_constr0(i)
532 theta_drange(i)=deg2rad*theta_drange(i)
534 C if(me.eq.king.or..not.out1file)
535 C & write (iout,*) 'FTORS',ftors
536 C do i=1,ntheta_constr
537 C ii = itheta_constr(i)
538 C thetabound(1,ii) = phi0(i)-drange(i)
539 C thetabound(2,ii) = phi0(i)+drange(i)
541 endif ! ntheta_constr.gt.0
542 endif! with_theta_constr
543 if (constr_homology.gt.0) then
544 c write (iout,*) "About to call read_constr_homology"
546 call read_constr_homology
547 c write (iout,*) "Exit read_constr_homology"
549 if (indpdb.gt.0 .or. pdbref) then
553 cref(j,i)=crefjlee(j,i)
558 write (iout,*) "Array C"
560 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
561 & (c(j,i+nres),j=1,3)
563 write (iout,*) "Array Cref"
565 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
566 & (cref(j,i+nres),j=1,3)
570 call int_from_cart1(.false.)
571 call sc_loc_geom(.false.)
575 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
579 dc(j,i)=c(j,i+1)-c(j,i)
580 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
585 dc(j,i+nres)=c(j,i+nres)-c(j,i)
586 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
593 write (iout,*) "calling read_saxs_consrtr",nsaxs
594 if (nsaxs.gt.0) call read_saxs_constr
597 c read(inp,'(a)') pdbfile
598 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
599 c open(ipdbin,file=pdbfile,status='old',err=33)
601 c 33 write (iout,'(a)') 'Error opening PDB file.'
604 c print *,'Begin reading pdb data'
606 c print *,'Finished reading pdb data'
607 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
609 c itype_pdb(i)=itype(i)
612 c write (iout,'(a,i3)') 'nsup=',nsup
614 c if (nsup.le.(nct-nnt+1)) then
615 c do i=0,nct-nnt+1-nsup
616 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
622 c & 'Error - sequences to be superposed do not match.'
625 c do i=0,nsup-(nct-nnt+1)
626 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
628 c nstart_sup=nstart_sup+i
634 c & 'Error - sequences to be superposed do not match.'
637 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
638 c & ' nstart_seq=',nstart_seq
643 C write (iout,'(/a,i3,a)')
644 C 'The chain contains',ns,' disulfide-bridging cysteines.'
645 write (iout,'(20i4)') (iss(i),i=1,ns)
647 write(iout,*)"Running with dynamic disulfide-bond formation"
649 write (iout,'(/a/)') 'Pre-formed links are:'
655 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
656 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
662 if (ns.gt.0.and.dyn_ss) then
666 forcon(i-nss)=forcon(i)
673 dyn_ss_mask(iss(i))=.true.
676 c Read distance restraints
677 if (constr_dist.gt.0) then
678 call read_dist_constr
683 c-----------------------------------------------------------------------------
684 logical function seq_comp(itypea,itypeb,length)
686 integer length,itypea(length),itypeb(length)
689 if (itypea(i).ne.itypeb(i)) then
697 c-----------------------------------------------------------------------------
698 subroutine read_bridge
699 C Read information about disulfide bridges.
702 include 'COMMON.IOUNITS'
705 include 'COMMON.INTERACT'
706 include 'COMMON.LOCAL'
707 include 'COMMON.NAMES'
708 include 'COMMON.CHAIN'
709 include 'COMMON.FFIELD'
710 include 'COMMON.SBRIDGE'
711 include 'COMMON.HEADER'
712 include 'COMMON.CONTROL'
713 include 'COMMON.TIME1'
715 include 'COMMON.INFO'
718 C Read bridging residues.
719 read (inp,*) ns,(iss(i),i=1,ns)
721 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
727 C Check whether the specified bridging residues are cystines.
729 if (itype(iss(i)).ne.1) then
730 write (iout,'(2a,i3,a)')
731 & 'Do you REALLY think that the residue ',
732 & restyp(itype(iss(i))),i,
733 & ' can form a disulfide bridge?!!!'
734 write (*,'(2a,i3,a)')
735 & 'Do you REALLY think that the residue ',
736 & restyp(itype(iss(i))),i,
737 & ' can form a disulfide bridge?!!!'
739 call mp_stopall(error_msg)
745 C Read preformed bridges.
747 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
750 C Check if the residues involved in bridges are in the specified list of
754 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
755 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
756 write (iout,'(a,i3,a)') 'Disulfide pair',i,
757 & ' contains residues present in other pairs.'
758 write (*,'(a,i3,a)') 'Disulfide pair',i,
759 & ' contains residues present in other pairs.'
761 call mp_stopall(error_msg)
768 if (ihpb(i).eq.iss(j)) goto 10
770 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
773 if (jhpb(i).eq.iss(j)) goto 20
775 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
788 c----------------------------------------------------------------------------
789 subroutine read_angles(kanal,*)
794 include 'COMMON.CHAIN'
795 include 'COMMON.IOUNITS'
797 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
798 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
799 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
800 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
802 theta(i)=deg2rad*theta(i)
803 phi(i)=deg2rad*phi(i)
804 alph(i)=deg2rad*alph(i)
805 omeg(i)=deg2rad*omeg(i)
810 c----------------------------------------------------------------------------
811 subroutine reada(rekord,lancuch,wartosc,default)
813 character*(*) rekord,lancuch
814 double precision wartosc,default
817 iread=index(rekord,lancuch)
822 iread=iread+ilen(lancuch)+1
823 read (rekord(iread:),*) wartosc
826 c----------------------------------------------------------------------------
827 subroutine multreada(rekord,lancuch,tablica,dim,default)
830 double precision tablica(dim),default
831 character*(*) rekord,lancuch
837 iread=index(rekord,lancuch)
838 if (iread.eq.0) return
839 iread=iread+ilen(lancuch)+1
840 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
843 c----------------------------------------------------------------------------
844 subroutine readi(rekord,lancuch,wartosc,default)
846 character*(*) rekord,lancuch
847 integer wartosc,default
850 iread=index(rekord,lancuch)
855 iread=iread+ilen(lancuch)+1
856 read (rekord(iread:),*) wartosc
859 C----------------------------------------------------------------------
860 subroutine multreadi(rekord,lancuch,tablica,dim,default)
863 integer tablica(dim),default
864 character*(*) rekord,lancuch
871 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
872 if (iread.eq.0) return
873 iread=iread+ilen(lancuch)+1
874 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
878 c----------------------------------------------------------------------------
879 subroutine card_concat(card)
881 include 'COMMON.IOUNITS'
883 character*80 karta,ucase
885 read (inp,'(a)') karta
888 do while (karta(80:80).eq.'&')
889 card=card(:ilen(card)+1)//karta(:79)
890 read (inp,'(a)') karta
893 card=card(:ilen(card)+1)//karta
896 c----------------------------------------------------------------------------
905 include 'COMMON.IOUNITS'
906 include 'COMMON.CONTROL'
907 integer lenpre,lenpot,ilen
909 character*16 cformat,cprint
911 integer lenint,lenout
912 call getenv('INPUT',prefix)
913 call getenv('OUTPUT',prefout)
914 call getenv('INTIN',prefintin)
915 call getenv('COORD',cformat)
916 call getenv('PRINTCOOR',cprint)
917 call getenv('SCRATCHDIR',scratchdir)
920 if (index(ucase(cformat),'CX').gt.0) then
927 lenint=ilen(prefintin)
928 C Get the names and open the input files
929 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
931 write (liczba,'(bz,i3.3)') me
932 outname=prefout(:lenout)//'_clust.out_'//liczba
934 outname=prefout(:lenout)//'_clust.out'
937 intinname=prefintin(:lenint)//'.bx'
938 else if (from_cx) then
939 intinname=prefintin(:lenint)//'.cx'
941 intinname=prefintin(:lenint)//'.int'
943 rmsname=prefintin(:lenint)//'.rms'
944 open (jplot,file=prefout(:ilen(prefout))//'.tex',
946 open (jrms,file=rmsname,status='unknown')
947 open(iout,file=outname,status='unknown')
948 C Get parameter filenames and open the parameter files.
949 call getenv('BONDPAR',bondname)
950 open (ibond,file=bondname,status='old')
951 call getenv('THETPAR',thetname)
952 open (ithep,file=thetname,status='old')
953 call getenv('ROTPAR',rotname)
954 open (irotam,file=rotname,status='old')
955 call getenv('TORPAR',torname)
956 open (itorp,file=torname,status='old')
958 call getenv('TORDPAR',tordname)
959 open (itordp,file=tordname,status='old')
961 call getenv('FOURIER',fouriername)
962 open (ifourier,file=fouriername,status='old')
963 call getenv('ELEPAR',elename)
964 open (ielep,file=elename,status='old')
965 call getenv('SIDEPAR',sidename)
966 open (isidep,file=sidename,status='old')
967 call getenv('SIDEP',sidepname)
968 open (isidep1,file=sidepname,status="old")
969 call getenv('SCCORPAR',sccorname)
970 open (isccor,file=sccorname,status="old")
971 call getenv('LIPTRANPAR',liptranname)
972 open (iliptranpar,file=liptranname,status='old')
975 C 8/9/01 In the newest version SCp interaction constants are read from a file
976 C Use -DOLDSCP to use hard-coded constants instead.
978 call getenv('SCPPAR',scpname)
979 open (iscpp,file=scpname,status='old')
983 c--------------------------------------------------------------------------
984 subroutine read_dist_constr
985 implicit real*8 (a-h,o-z)
990 include 'COMMON.CONTROL'
991 include 'COMMON.CHAIN'
992 include 'COMMON.IOUNITS'
993 include 'COMMON.SBRIDGE'
994 include 'COMMON.INTERACT'
995 integer ifrag_(2,100),ipair_(2,100)
996 double precision wfrag_(100),wpair_(100)
997 character*500 controlcard
998 logical lprn /.true./
999 logical normalize,next
1001 double precision scal_bfac
1002 double precision xlink(4,0:4) /
1004 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
1005 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
1006 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
1007 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
1008 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
1009 write (iout,*) "Calling read_dist_constr"
1010 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1016 call card_concat(controlcard)
1017 next = index(controlcard,"NEXT").gt.0
1018 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
1019 write (iout,*) "restr_type",restr_type
1020 call readi(controlcard,"NFRAG",nfrag_,0)
1021 call readi(controlcard,"NFRAG",nfrag_,0)
1022 call readi(controlcard,"NPAIR",npair_,0)
1023 call readi(controlcard,"NDIST",ndist_,0)
1024 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1025 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
1026 if (restr_type.eq.10)
1027 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
1028 if (restr_type.eq.12)
1029 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
1030 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1031 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1032 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1033 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1034 normalize = index(controlcard,"NORMALIZE").gt.0
1035 write (iout,*) "WBOLTZD",wboltzd
1036 write (iout,*) "SCAL_PEAK",scal_peak
1037 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1038 write (iout,*) "IFRAG"
1040 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1042 write (iout,*) "IPAIR"
1044 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1046 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
1048 read(inp,'(a)') pdbfile
1050 & "Distance restraints will be constructed from structure ",pdbfile
1051 open(ipdbin,file=pdbfile,status='old',err=11)
1052 call readpdb(.true.)
1057 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1058 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1059 & ifrag_(2,i)=nstart_sup+nsup-1
1060 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1062 if (wfrag_(i).eq.0.0d0) cycle
1063 do j=ifrag_(1,i),ifrag_(2,i)-1
1064 do k=j+1,ifrag_(2,i)
1065 c write (iout,*) "j",j," k",k
1067 if (restr_type.eq.1) then
1073 forcon(nhpb)=wfrag_(i)
1074 else if (constr_dist.eq.2) then
1075 if (ddjk.le.dist_cut) then
1081 forcon(nhpb)=wfrag_(i)
1083 else if (restr_type.eq.3) then
1089 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1091 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
1092 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1097 if (wpair_(i).eq.0.0d0) cycle
1105 do j=ifrag_(1,ii),ifrag_(2,ii)
1106 do k=ifrag_(1,jj),ifrag_(2,jj)
1108 if (restr_type.eq.1) then
1114 forcon(nhpb)=wpair_(i)
1115 else if (constr_dist.eq.2) then
1116 if (ddjk.le.dist_cut) then
1122 forcon(nhpb)=wpair_(i)
1124 else if (restr_type.eq.3) then
1130 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1132 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1133 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1139 write (iout,*) "Distance restraints as read from input"
1141 if (restr_type.eq.12) then
1142 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1143 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1144 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1145 & fordepth_peak(nhpb_peak+1),npeak
1146 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1147 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1148 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1149 c & fordepth_peak(nhpb_peak+1),npeak
1150 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1151 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1152 nhpb_peak=nhpb_peak+1
1153 irestr_type_peak(nhpb_peak)=12
1154 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1156 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1157 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1158 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1159 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1160 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1161 if (ibecarb_peak(nhpb_peak).eq.3) then
1162 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1163 else if (ibecarb_peak(nhpb_peak).eq.2) then
1164 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1165 else if (ibecarb_peak(nhpb_peak).eq.1) then
1166 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1167 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1169 else if (restr_type.eq.11) then
1170 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1171 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1172 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1173 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1175 irestr_type(nhpb)=11
1176 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1177 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1178 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1179 c if (ibecarb(nhpb).gt.0) then
1180 c ihpb(nhpb)=ihpb(nhpb)+nres
1181 c jhpb(nhpb)=jhpb(nhpb)+nres
1183 if (ibecarb(nhpb).eq.3) then
1184 ihpb(nhpb)=ihpb(nhpb)+nres
1185 else if (ibecarb(nhpb).eq.2) then
1186 ihpb(nhpb)=ihpb(nhpb)+nres
1187 else if (ibecarb(nhpb).eq.1) then
1188 ihpb(nhpb)=ihpb(nhpb)+nres
1189 jhpb(nhpb)=jhpb(nhpb)+nres
1191 else if (restr_type.eq.10) then
1192 c Cross-lonk Markov-like potential
1193 call card_concat(controlcard)
1194 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1195 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1197 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1198 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1199 if (index(controlcard,"ZL").gt.0) then
1201 else if (index(controlcard,"ADH").gt.0) then
1203 else if (index(controlcard,"PDH").gt.0) then
1205 else if (index(controlcard,"DSS").gt.0) then
1210 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1211 & xlink(1,link_type))
1212 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1213 & xlink(2,link_type))
1214 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1215 & xlink(3,link_type))
1216 call reada(controlcard,"SIGMA",forcon(nhpb+1),
1217 & xlink(4,link_type))
1218 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1219 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1220 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1221 if (forcon(nhpb+1).le.0.0d0 .or.
1222 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1224 irestr_type(nhpb)=10
1225 c if (ibecarb(nhpb).gt.0) then
1226 c ihpb(nhpb)=ihpb(nhpb)+nres
1227 c jhpb(nhpb)=jhpb(nhpb)+nres
1229 if (ibecarb(nhpb).eq.3) then
1230 jhpb(nhpb)=jhpb(nhpb)+nres
1231 else if (ibecarb(nhpb).eq.2) then
1232 ihpb(nhpb)=ihpb(nhpb)+nres
1233 else if (ibecarb(nhpb).eq.1) then
1234 ihpb(nhpb)=ihpb(nhpb)+nres
1235 jhpb(nhpb)=jhpb(nhpb)+nres
1237 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1238 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1239 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1243 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1244 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1245 if (forcon(nhpb+1).gt.0.0d0) then
1247 if (dhpb1(nhpb).eq.0.0d0) then
1252 c if (ibecarb(nhpb).gt.0) then
1253 c ihpb(nhpb)=ihpb(nhpb)+nres
1254 c jhpb(nhpb)=jhpb(nhpb)+nres
1256 if (ibecarb(nhpb).eq.3) then
1257 jhpb(nhpb)=jhpb(nhpb)+nres
1258 else if (ibecarb(nhpb).eq.2) then
1259 ihpb(nhpb)=ihpb(nhpb)+nres
1260 else if (ibecarb(nhpb).eq.1) then
1261 ihpb(nhpb)=ihpb(nhpb)+nres
1262 jhpb(nhpb)=jhpb(nhpb)+nres
1264 if (dhpb(nhpb).eq.0.0d0)
1265 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1267 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1268 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1270 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1271 C if (forcon(nhpb+1).gt.0.0d0) then
1273 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1276 if (restr_type.eq.4) then
1277 write (iout,*) "The BFAC array"
1279 write (iout,'(i5,f10.5)') i,bfac(i)
1282 if (itype(i).eq.ntyp1) cycle
1284 if (itype(j).eq.ntyp1) cycle
1285 if (itype(i).eq.10) then
1290 if (itype(j).eq.10) then
1300 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1303 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1304 ihpb(nhpb)=i+nres*ii
1305 jhpb(nhpb)=j+nres*jj
1306 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1307 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1308 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1309 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1318 if (restr_type.eq.5) then
1319 restr_on_coord=.true.
1321 if (itype(i).eq.ntyp1) cycle
1322 bfac(i)=(scal_bfac/bfac(i))**2
1331 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1332 & fordepthmax=fordepth(i)
1335 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1338 if (nhpb.gt.nss) then
1339 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1340 & "The following",nhpb-nss,
1341 & " distance restraints have been imposed:",
1342 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1345 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1346 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1353 11 write (iout,*)"read_dist_restr: error reading reference structure"
1356 c-------------------------------------------------------------------------------
1357 subroutine read_saxs_constr
1358 implicit real*8 (a-h,o-z)
1359 include 'DIMENSIONS'
1363 include 'COMMON.CONTROL'
1364 include 'COMMON.CHAIN'
1365 include 'COMMON.IOUNITS'
1366 include 'COMMON.SBRIDGE'
1367 include 'COMMON.SAXS'
1368 double precision cm(3)
1370 write (iout,*) "Calling read_saxs nsaxs",nsaxs
1372 if (saxs_mode.eq.0) then
1373 c SAXS distance distribution
1375 read(inp,*) distsaxs(i),Psaxs(i)
1379 Cnorm = Cnorm + Psaxs(i)
1381 write (iout,*) "Cnorm",Cnorm
1383 Psaxs(i)=Psaxs(i)/Cnorm
1385 write (iout,*) "Normalized distance distribution from SAXS"
1387 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1391 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1393 write (iout,*) "Wsaxs0",Wsaxs0
1397 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1404 cm(j)=cm(j)+Csaxs(j,i)
1412 Csaxs(j,i)=Csaxs(j,i)-cm(j)
1415 write (iout,*) "SAXS sphere coordinates"
1417 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)