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 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
31 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
32 call readi(controlcard,'TORMODE',tor_mode,0)
33 call readi(controlcard,'NRES',nres,0)
34 call readi(controlcard,'RESCALE',rescale_mode,2)
35 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
36 write (iout,*) "DISTCHAINMAX",distchainmax
37 C Reading the dimensions of box in x,y,z coordinates
38 call reada(controlcard,'BOXX',boxxsize,100.0d0)
39 call reada(controlcard,'BOXY',boxysize,100.0d0)
40 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
41 c Cutoff range for interactions
42 call reada(controlcard,"R_CUT",r_cut,25.0d0)
43 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
44 write (iout,*) "Cutoff on interactions",r_cut
45 write (iout,*) "lambda",rlamb
46 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
47 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
48 if (lipthick.gt.0.0d0) then
49 bordliptop=(boxzsize+lipthick)/2.0
50 bordlipbot=bordliptop-lipthick
52 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
53 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
54 buflipbot=bordlipbot+lipbufthick
55 bufliptop=bordliptop-lipbufthick
56 if ((lipbufthick*2.0d0).gt.lipthick)
57 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
59 write(iout,*) "bordliptop=",bordliptop
60 write(iout,*) "bordlipbot=",bordlipbot
61 write(iout,*) "bufliptop=",bufliptop
62 write(iout,*) "buflipbot=",buflipbot
64 call readi(controlcard,'SHIELD',shield_mode,0)
65 write (iout,*) "SHIELD MODE",shield_mode
66 if (shield_mode.gt.0) then
67 pdbref=(index(controlcard,'PDBREF').gt.0)
68 if (index(controlcard,"CASC").gt.0) then
70 else if (index(controlcard,"SCONLY").gt.0) then
76 C VSolvSphere the volume of solving sphere
78 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
79 C there will be no distinction between proline peptide group and normal peptide
80 C group in case of shielding parameters
81 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
82 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
83 write (iout,*) VSolvSphere,VSolvSphere_div
84 C long axis of side chain
86 long_r_sidechain(i)=vbldsc0(1,i)
87 short_r_sidechain(i)=sigma0(i)
91 call readi(controlcard,'PDBOUT',outpdb,0)
92 call readi(controlcard,'MOL2OUT',outmol2,0)
93 refstr=(index(controlcard,'REFSTR').gt.0)
94 pdbref=(index(controlcard,'PDBREF').gt.0)
95 refstr = refstr .or. pdbref
96 write (iout,*) "REFSTR",refstr," PDBREF",pdbref
97 iscode=index(controlcard,'ONE_LETTER')
98 tree=(index(controlcard,'MAKE_TREE').gt.0)
99 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 write (iout,*) "with_dihed_constr ",with_dihed_constr,
102 & " CONSTR_DIST",constr_dist
103 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
104 write (iout,*) "with_theta_constr ",with_theta_constr
106 min_var=(index(controlcard,'MINVAR').gt.0)
107 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
108 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
109 print_fittest=(index(controlcard,'PRINT_FITTEST').gt.0)
110 call readi(controlcard,'NCUT',ncut,0)
112 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
115 call readi(controlcard,'NCLUST',nclust,5)
117 call readi(controlcard,'SYM',symetr,1)
118 write (iout,*) 'sym', symetr
119 call readi(controlcard,'NSTART',nstart,0)
120 call readi(controlcard,'NEND',nend,0)
121 call reada(controlcard,'ECUT',ecut,10.0d0)
122 call reada(controlcard,'PROB',prob_limit,0.99d0)
123 write (iout,*) "Probability limit",prob_limit
124 lgrp=(index(controlcard,'LGRP').gt.0)
125 caonly=(index(controlcard,'CA_ONLY').gt.0)
126 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
127 call readi(controlcard,'IOPT',iopt,2)
128 lefree = index(controlcard,"EFREE").gt.0
129 call readi(controlcard,'NTEMP',nT,1)
130 write (iout,*) "nT",nT
131 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
132 write (iout,*) "nT",nT
133 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
135 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
137 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
138 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
139 lprint_int=index(controlcard,"PRINT_INT") .gt.0
140 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
141 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
142 c if (constr_homology) tole=dmax1(tole,1.5d0)
143 write (iout,*) "with_homology_constr ",with_dihed_constr,
144 & " CONSTR_HOMOLOGY",constr_homology
145 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
146 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
147 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
148 write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD
149 call readi(controlcard,'NSAXS',nsaxs,0)
150 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
151 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
152 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
153 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
154 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
158 c--------------------------------------------------------------------------
161 C Read molecular data.
165 include 'COMMON.IOUNITS'
168 include 'COMMON.INTERACT'
169 include 'COMMON.LOCAL'
170 include 'COMMON.NAMES'
171 include 'COMMON.CHAIN'
172 include 'COMMON.FFIELD'
173 include 'COMMON.SBRIDGE'
174 include 'COMMON.HEADER'
175 include 'COMMON.CONTROL'
176 include 'COMMON.CONTACTS'
177 include 'COMMON.TIME1'
178 include 'COMMON.TORCNSTR'
179 include 'COMMON.SHIELD'
180 include 'COMMON.SAXS'
182 include 'COMMON.INFO'
184 character*4 sequence(maxres)
185 character*800 weightcard,controlcard
187 double precision x(maxvar)
188 double precision phihel,phibet,sigmahel,sigmabet,sumv,
190 integer itype_pdb(maxres)
192 integer i,j,kkk,i1,i2,it1,it2,tperm,ii,iperm
196 C Read weights of the subsequent energy terms.
197 call card_concat(weightcard)
198 call reada(weightcard,'WSC',wsc,1.0d0)
199 call reada(weightcard,'WLONG',wsc,wsc)
200 call reada(weightcard,'WSCP',wscp,1.0d0)
201 call reada(weightcard,'WELEC',welec,1.0D0)
202 call reada(weightcard,'WVDWPP',wvdwpp,welec)
203 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
204 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
205 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
206 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
207 call reada(weightcard,'WTURN3',wturn3,1.0D0)
208 call reada(weightcard,'WTURN4',wturn4,1.0D0)
209 call reada(weightcard,'WTURN6',wturn6,1.0D0)
210 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
211 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
212 call reada(weightcard,'WBOND',wbond,1.0D0)
213 call reada(weightcard,'WTOR',wtor,1.0D0)
214 call reada(weightcard,'WTORD',wtor_d,1.0D0)
215 call reada(weightcard,'WANG',wang,1.0D0)
216 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
217 call reada(weightcard,'WSAXS',wsaxs,0.0D0)
218 call reada(weightcard,'SCAL14',scal14,0.4D0)
219 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
220 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
221 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
222 if (index(weightcard,'SOFT').gt.0) ipot=6
223 call reada(weightcard,"D0CM",d0cm,3.78d0)
224 call reada(weightcard,"AKCM",akcm,15.1d0)
225 call reada(weightcard,"AKTH",akth,11.0d0)
226 call reada(weightcard,"AKCT",akct,12.0d0)
227 call reada(weightcard,"V1SS",v1ss,-1.08d0)
228 call reada(weightcard,"V2SS",v2ss,7.61d0)
229 call reada(weightcard,"V3SS",v3ss,13.7d0)
230 call reada(weightcard,"EBR",ebr,-5.50D0)
231 call reada(weightcard,'WSHIELD',wshield,1.0d0)
232 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
233 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
234 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
235 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
236 call reada(weightcard,'WLT',wliptran,0.0D0)
237 call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
238 call reada(weightcard,"ATRISS",atriss,0.301D0)
239 call reada(weightcard,"BTRISS",btriss,0.021D0)
240 call reada(weightcard,"CTRISS",ctriss,1.001D0)
241 call reada(weightcard,"DTRISS",dtriss,1.001D0)
242 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
244 dyn_ss_mask(i)=.false.
248 dyn_ssbond_ij(i,j)=1.0d300
251 call reada(weightcard,"HT",Ht,0.0D0)
253 ss_depth=ebr/wsc-0.25*eps(1,1)
254 Ht=Ht/wsc-0.25*eps(1,1)
255 akcm=akcm*wstrain/wsc
256 akth=akth*wstrain/wsc
257 akct=akct*wstrain/wsc
258 v1ss=v1ss*wstrain/wsc
259 v2ss=v2ss*wstrain/wsc
260 v3ss=v3ss*wstrain/wsc
262 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
264 write (iout,'(/a)') "Disulfide bridge parameters:"
265 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
266 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
267 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
268 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
270 write (iout,*) "Parameters of the 'trisulfide' potential"
271 write (iout,*) "ATRISS=", atriss
272 write (iout,*) "BTRISS=", btriss
273 write (iout,*) "CTRISS=", ctriss
274 write (iout,*) "DTRISS=", dtriss
276 C 12/1/95 Added weight for the multi-body term WCORR
277 call reada(weightcard,'WCORRH',wcorr,1.0D0)
278 if (wcorr4.gt.0.0d0) wcorr=wcorr4
298 weights(28)=wdfa_dist
301 weights(31)=wdfa_beta
302 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
303 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
304 & wturn4,wturn6,wsccor
305 10 format (/'Energy-term weights (unscaled):'//
306 & 'WSCC= ',f10.6,' (SC-SC)'/
307 & 'WSCP= ',f10.6,' (SC-p)'/
308 & 'WELEC= ',f10.6,' (p-p electr)'/
309 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
310 & 'WBOND= ',f10.6,' (stretching)'/
311 & 'WANG= ',f10.6,' (bending)'/
312 & 'WSCLOC= ',f10.6,' (SC local)'/
313 & 'WTOR= ',f10.6,' (torsional)'/
314 & 'WTORD= ',f10.6,' (double torsional)'/
315 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
316 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
317 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
318 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
319 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
320 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
321 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
322 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
323 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
325 if (wcorr4.gt.0.0d0) then
326 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
327 & 'between contact pairs of peptide groups'
328 write (iout,'(2(a,f5.3/))')
329 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
330 & 'Range of quenching the correlation terms:',2*delt_corr
331 else if (wcorr.gt.0.0d0) then
332 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
333 & 'between contact pairs of peptide groups'
335 write (iout,'(a,f8.3)')
336 & 'Scaling factor of 1,4 SC-p interactions:',scal14
337 write (iout,'(a,f8.3)')
338 & 'General scaling factor of SC-p interactions:',scalscp
339 r0_corr=cutoff_corr-delt_corr
341 aad(i,1)=scalscp*aad(i,1)
342 aad(i,2)=scalscp*aad(i,2)
343 bad(i,1)=scalscp*bad(i,1)
344 bad(i,2)=scalscp*bad(i,2)
347 write (iout,'(/a/)') "DFA pseudopotential parameters:"
348 write (iout,'(a,f10.6,a)')
349 & "WDFAD= ",wdfa_dist," (distance)",
350 & "WDFAT= ",wdfa_tor," (backbone angles)",
351 & "WDFAN= ",wdfa_nei," (neighbors)",
352 & "WDFAB= ",wdfa_beta," (beta structure)"
355 c print *,'indpdb=',indpdb,' pdbref=',pdbref
357 C Read sequence if not taken from the pdb file.
358 if (iscode.gt.0) then
359 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
361 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
363 C Convert sequence to numeric code
365 itype(i)=rescode(i,sequence(i),iscode)
368 c print '(20i4)',(itype(i),i=1,nres)
372 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
374 if (itype(i).eq.ntyp1) then
378 else if (iabs(itype(i+1)).ne.20) then
380 else if (iabs(itype(i)).ne.20) then
387 write (iout,*) "ITEL"
389 write (iout,*) i,itype(i),itel(i)
392 c print *,'Call Read_Bridge.'
394 C this fragment reads diheadral constrains
397 c print *,'NNT=',NNT,' NCT=',NCT
398 call seq2chains(nres,itype,nchain,chain_length,chain_border,
400 write(iout,*) "nres",nres," nchain",nchain
402 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
405 call chain_symmetry(nchain,nres,itype,chain_border,
406 & chain_length,npermchain,tabpermchain)
408 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
411 write(iout,*) "residue permutations"
413 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
415 if (itype(1).eq.ntyp1) nnt=2
416 if (itype(nres).eq.ntyp1) nct=nct-1
417 if (nstart.lt.nnt) nstart=nnt
418 if (nend.gt.nct .or. nend.eq.0) nend=nct
419 write (iout,*) "nstart",nstart," nend",nend
422 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
423 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
425 print*, 'init_dfa_vars finished!'
427 print*, 'read_dfa_info finished!'
430 C If the reference structure is not read set the superposition
437 if (with_dihed_constr) then
439 read (inp,*) ndih_constr
440 if (ndih_constr.gt.0) then
443 C write (iout,*) 'FTORS',ftors
444 C ftors is the force constant for torsional quartic constrains
445 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
448 & 'There are',ndih_constr,' constraints on phi angles.'
450 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
454 phi0(i)=deg2rad*phi0(i)
455 drange(i)=deg2rad*drange(i)
457 else if (ndih_constr.lt.0) then
459 call card_concat(controlcard)
460 call reada(controlcard,"PHIHEL",phihel,50.0D0)
461 call reada(controlcard,"PHIBET",phibet,180.0D0)
462 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
463 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
464 call reada(controlcard,"WDIHC",wdihc,0.591d0)
465 write (iout,*) "Weight of the dihedral restraint term",wdihc
466 read(inp,'(9x,3f7.3)')
467 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
468 write (iout,*) "The secprob array"
470 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
474 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
475 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
476 ndih_constr=ndih_constr+1
477 idih_constr(ndih_constr)=i
480 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
481 sumv=sumv+vpsipred(j,ndih_constr)
484 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
486 phibound(1,ndih_constr)=phihel*deg2rad
487 phibound(2,ndih_constr)=phibet*deg2rad
488 sdihed(1,ndih_constr)=sigmahel*deg2rad
489 sdihed(2,ndih_constr)=sigmabet*deg2rad
493 & 'There are',ndih_constr,
494 & ' bimodal restraints on gamma angles.'
496 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
497 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
498 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
499 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
500 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
501 & (vpsipred(j,i),j=1,3)
504 endif ! endif ndif_constr.gt.0
505 endif ! with_dihed_constr
506 if (with_theta_constr) then
507 C with_theta_constr is keyword allowing for occurance of theta constrains
508 read (inp,*) ntheta_constr
509 C ntheta_constr is the number of theta constrains
510 if (ntheta_constr.gt.0) then
512 read (inp,*) (itheta_constr(i),theta_constr0(i),
513 & theta_drange(i),for_thet_constr(i),
515 C the above code reads from 1 to ntheta_constr
516 C itheta_constr(i) residue i for which is theta_constr
517 C theta_constr0 the global minimum value
518 C theta_drange is range for which there is no energy penalty
519 C for_thet_constr is the force constant for quartic energy penalty
522 & 'There are',ntheta_constr,' constraints on phi angles.'
524 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
530 theta_constr0(i)=deg2rad*theta_constr0(i)
531 theta_drange(i)=deg2rad*theta_drange(i)
533 C if(me.eq.king.or..not.out1file)
534 C & write (iout,*) 'FTORS',ftors
535 C do i=1,ntheta_constr
536 C ii = itheta_constr(i)
537 C thetabound(1,ii) = phi0(i)-drange(i)
538 C thetabound(2,ii) = phi0(i)+drange(i)
540 endif ! ntheta_constr.gt.0
541 endif! with_theta_constr
542 if (constr_homology.gt.0) then
543 c write (iout,*) "About to call read_constr_homology"
545 call read_constr_homology
546 c write (iout,*) "Exit read_constr_homology"
548 if (indpdb.gt.0 .or. pdbref) then
552 cref(j,i)=crefjlee(j,i)
557 write (iout,*) "Array C"
559 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
560 & (c(j,i+nres),j=1,3)
562 write (iout,*) "Array Cref"
564 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
565 & (cref(j,i+nres),j=1,3)
569 call int_from_cart1(.false.)
570 call sc_loc_geom(.false.)
574 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
578 dc(j,i)=c(j,i+1)-c(j,i)
579 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
584 dc(j,i+nres)=c(j,i+nres)-c(j,i)
585 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
592 write (iout,*) "calling read_saxs_consrtr",nsaxs
593 if (nsaxs.gt.0) call read_saxs_constr
596 c read(inp,'(a)') pdbfile
597 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
598 c open(ipdbin,file=pdbfile,status='old',err=33)
600 c 33 write (iout,'(a)') 'Error opening PDB file.'
603 c print *,'Begin reading pdb data'
605 c print *,'Finished reading pdb data'
606 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
608 c itype_pdb(i)=itype(i)
611 c write (iout,'(a,i3)') 'nsup=',nsup
613 c if (nsup.le.(nct-nnt+1)) then
614 c do i=0,nct-nnt+1-nsup
615 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
621 c & 'Error - sequences to be superposed do not match.'
624 c do i=0,nsup-(nct-nnt+1)
625 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
627 c nstart_sup=nstart_sup+i
633 c & 'Error - sequences to be superposed do not match.'
636 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
637 c & ' nstart_seq=',nstart_seq
642 C write (iout,'(/a,i3,a)')
643 C 'The chain contains',ns,' disulfide-bridging cysteines.'
644 write (iout,'(20i4)') (iss(i),i=1,ns)
646 write(iout,*)"Running with dynamic disulfide-bond formation"
648 write (iout,'(/a/)') 'Pre-formed links are:'
654 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
655 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
661 if (ns.gt.0.and.dyn_ss) then
665 forcon(i-nss)=forcon(i)
672 dyn_ss_mask(iss(i))=.true.
675 c Read distance restraints
676 if (constr_dist.gt.0) then
677 call read_dist_constr
682 c-----------------------------------------------------------------------------
683 logical function seq_comp(itypea,itypeb,length)
685 integer length,itypea(length),itypeb(length)
688 if (itypea(i).ne.itypeb(i)) then
696 c-----------------------------------------------------------------------------
697 subroutine read_bridge
698 C Read information about disulfide bridges.
701 include 'COMMON.IOUNITS'
704 include 'COMMON.INTERACT'
705 include 'COMMON.LOCAL'
706 include 'COMMON.NAMES'
707 include 'COMMON.CHAIN'
708 include 'COMMON.FFIELD'
709 include 'COMMON.SBRIDGE'
710 include 'COMMON.HEADER'
711 include 'COMMON.CONTROL'
712 include 'COMMON.TIME1'
714 include 'COMMON.INFO'
717 C Read bridging residues.
718 read (inp,*) ns,(iss(i),i=1,ns)
720 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
726 C Check whether the specified bridging residues are cystines.
728 if (itype(iss(i)).ne.1) then
729 write (iout,'(2a,i3,a)')
730 & 'Do you REALLY think that the residue ',
731 & restyp(itype(iss(i))),i,
732 & ' can form a disulfide bridge?!!!'
733 write (*,'(2a,i3,a)')
734 & 'Do you REALLY think that the residue ',
735 & restyp(itype(iss(i))),i,
736 & ' can form a disulfide bridge?!!!'
738 call mp_stopall(error_msg)
744 C Read preformed bridges.
746 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
749 C Check if the residues involved in bridges are in the specified list of
753 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
754 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
755 write (iout,'(a,i3,a)') 'Disulfide pair',i,
756 & ' contains residues present in other pairs.'
757 write (*,'(a,i3,a)') 'Disulfide pair',i,
758 & ' contains residues present in other pairs.'
760 call mp_stopall(error_msg)
767 if (ihpb(i).eq.iss(j)) goto 10
769 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
772 if (jhpb(i).eq.iss(j)) goto 20
774 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
787 c----------------------------------------------------------------------------
788 subroutine read_angles(kanal,*)
793 include 'COMMON.CHAIN'
794 include 'COMMON.IOUNITS'
796 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
797 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
798 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
799 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
801 theta(i)=deg2rad*theta(i)
802 phi(i)=deg2rad*phi(i)
803 alph(i)=deg2rad*alph(i)
804 omeg(i)=deg2rad*omeg(i)
809 c----------------------------------------------------------------------------
810 subroutine reada(rekord,lancuch,wartosc,default)
812 character*(*) rekord,lancuch
813 double precision wartosc,default
816 iread=index(rekord,lancuch)
821 iread=iread+ilen(lancuch)+1
822 read (rekord(iread:),*) wartosc
825 c----------------------------------------------------------------------------
826 subroutine multreada(rekord,lancuch,tablica,dim,default)
829 double precision tablica(dim),default
830 character*(*) rekord,lancuch
836 iread=index(rekord,lancuch)
837 if (iread.eq.0) return
838 iread=iread+ilen(lancuch)+1
839 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
842 c----------------------------------------------------------------------------
843 subroutine readi(rekord,lancuch,wartosc,default)
845 character*(*) rekord,lancuch
846 integer wartosc,default
849 iread=index(rekord,lancuch)
854 iread=iread+ilen(lancuch)+1
855 read (rekord(iread:),*) wartosc
858 C----------------------------------------------------------------------
859 subroutine multreadi(rekord,lancuch,tablica,dim,default)
862 integer tablica(dim),default
863 character*(*) rekord,lancuch
870 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
871 if (iread.eq.0) return
872 iread=iread+ilen(lancuch)+1
873 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
877 c----------------------------------------------------------------------------
878 subroutine card_concat(card)
880 include 'COMMON.IOUNITS'
882 character*80 karta,ucase
884 read (inp,'(a)') karta
887 do while (karta(80:80).eq.'&')
888 card=card(:ilen(card)+1)//karta(:79)
889 read (inp,'(a)') karta
892 card=card(:ilen(card)+1)//karta
895 c----------------------------------------------------------------------------
904 include 'COMMON.IOUNITS'
905 include 'COMMON.CONTROL'
906 integer lenpre,lenpot,ilen
908 character*16 cformat,cprint
910 integer lenint,lenout
911 call getenv('INPUT',prefix)
912 call getenv('OUTPUT',prefout)
913 call getenv('INTIN',prefintin)
914 call getenv('COORD',cformat)
915 call getenv('PRINTCOOR',cprint)
916 call getenv('SCRATCHDIR',scratchdir)
919 if (index(ucase(cformat),'CX').gt.0) then
926 lenint=ilen(prefintin)
927 C Get the names and open the input files
928 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
930 write (liczba,'(bz,i3.3)') me
931 outname=prefout(:lenout)//'_clust.out_'//liczba
933 outname=prefout(:lenout)//'_clust.out'
936 intinname=prefintin(:lenint)//'.bx'
937 else if (from_cx) then
938 intinname=prefintin(:lenint)//'.cx'
940 intinname=prefintin(:lenint)//'.int'
942 rmsname=prefintin(:lenint)//'.rms'
943 open (jplot,file=prefout(:ilen(prefout))//'.tex',
945 open (jrms,file=rmsname,status='unknown')
946 open(iout,file=outname,status='unknown')
947 C Get parameter filenames and open the parameter files.
948 call getenv('BONDPAR',bondname)
949 open (ibond,file=bondname,status='old')
950 call getenv('THETPAR',thetname)
951 open (ithep,file=thetname,status='old')
952 call getenv('ROTPAR',rotname)
953 open (irotam,file=rotname,status='old')
954 call getenv('TORPAR',torname)
955 open (itorp,file=torname,status='old')
957 call getenv('TORDPAR',tordname)
958 open (itordp,file=tordname,status='old')
960 call getenv('FOURIER',fouriername)
961 open (ifourier,file=fouriername,status='old')
962 call getenv('ELEPAR',elename)
963 open (ielep,file=elename,status='old')
964 call getenv('SIDEPAR',sidename)
965 open (isidep,file=sidename,status='old')
966 call getenv('SIDEP',sidepname)
967 open (isidep1,file=sidepname,status="old")
968 call getenv('SCCORPAR',sccorname)
969 open (isccor,file=sccorname,status="old")
970 call getenv('LIPTRANPAR',liptranname)
971 open (iliptranpar,file=liptranname,status='old')
974 C 8/9/01 In the newest version SCp interaction constants are read from a file
975 C Use -DOLDSCP to use hard-coded constants instead.
977 call getenv('SCPPAR',scpname)
978 open (iscpp,file=scpname,status='old')
982 c--------------------------------------------------------------------------
983 subroutine read_dist_constr
984 implicit real*8 (a-h,o-z)
989 include 'COMMON.CONTROL'
990 include 'COMMON.CHAIN'
991 include 'COMMON.IOUNITS'
992 include 'COMMON.SBRIDGE'
993 include 'COMMON.INTERACT'
994 integer ifrag_(2,100),ipair_(2,100)
995 double precision wfrag_(100),wpair_(100)
996 character*500 controlcard
997 logical lprn /.true./
998 logical normalize,next
1000 double precision scal_bfac
1001 double precision xlink(4,0:4) /
1003 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
1004 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
1005 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
1006 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
1007 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
1008 write (iout,*) "Calling read_dist_constr"
1009 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1015 call card_concat(controlcard)
1016 next = index(controlcard,"NEXT").gt.0
1017 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
1018 write (iout,*) "restr_type",restr_type
1019 call readi(controlcard,"NFRAG",nfrag_,0)
1020 call readi(controlcard,"NFRAG",nfrag_,0)
1021 call readi(controlcard,"NPAIR",npair_,0)
1022 call readi(controlcard,"NDIST",ndist_,0)
1023 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1024 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
1025 if (restr_type.eq.10)
1026 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
1027 if (restr_type.eq.12)
1028 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
1029 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1030 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1031 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1032 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1033 normalize = index(controlcard,"NORMALIZE").gt.0
1034 write (iout,*) "WBOLTZD",wboltzd
1035 write (iout,*) "SCAL_PEAK",scal_peak
1036 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1037 write (iout,*) "IFRAG"
1039 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1041 write (iout,*) "IPAIR"
1043 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1045 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
1047 read(inp,'(a)') pdbfile
1049 & "Distance restraints will be constructed from structure ",pdbfile
1050 open(ipdbin,file=pdbfile,status='old',err=11)
1051 call readpdb(.true.)
1056 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1057 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1058 & ifrag_(2,i)=nstart_sup+nsup-1
1059 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1061 if (wfrag_(i).eq.0.0d0) cycle
1062 do j=ifrag_(1,i),ifrag_(2,i)-1
1063 do k=j+1,ifrag_(2,i)
1064 c write (iout,*) "j",j," k",k
1066 if (restr_type.eq.1) then
1072 forcon(nhpb)=wfrag_(i)
1073 else if (constr_dist.eq.2) then
1074 if (ddjk.le.dist_cut) then
1080 forcon(nhpb)=wfrag_(i)
1082 else if (restr_type.eq.3) then
1088 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1090 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
1091 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1096 if (wpair_(i).eq.0.0d0) cycle
1104 do j=ifrag_(1,ii),ifrag_(2,ii)
1105 do k=ifrag_(1,jj),ifrag_(2,jj)
1107 if (restr_type.eq.1) then
1113 forcon(nhpb)=wpair_(i)
1114 else if (constr_dist.eq.2) then
1115 if (ddjk.le.dist_cut) then
1121 forcon(nhpb)=wpair_(i)
1123 else if (restr_type.eq.3) then
1129 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1131 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1132 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1138 write (iout,*) "Distance restraints as read from input"
1140 if (restr_type.eq.12) then
1141 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1142 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1143 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1144 & fordepth_peak(nhpb_peak+1),npeak
1145 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1146 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1147 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1148 c & fordepth_peak(nhpb_peak+1),npeak
1149 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1150 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1151 nhpb_peak=nhpb_peak+1
1152 irestr_type_peak(nhpb_peak)=12
1153 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1155 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1156 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1157 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1158 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1159 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1160 if (ibecarb_peak(nhpb_peak).eq.3) then
1161 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1162 else if (ibecarb_peak(nhpb_peak).eq.2) then
1163 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1164 else if (ibecarb_peak(nhpb_peak).eq.1) then
1165 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1166 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1168 else if (restr_type.eq.11) then
1169 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1170 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1171 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1172 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1174 irestr_type(nhpb)=11
1175 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1176 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1177 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1178 c if (ibecarb(nhpb).gt.0) then
1179 c ihpb(nhpb)=ihpb(nhpb)+nres
1180 c jhpb(nhpb)=jhpb(nhpb)+nres
1182 if (ibecarb(nhpb).eq.3) then
1183 ihpb(nhpb)=ihpb(nhpb)+nres
1184 else if (ibecarb(nhpb).eq.2) then
1185 ihpb(nhpb)=ihpb(nhpb)+nres
1186 else if (ibecarb(nhpb).eq.1) then
1187 ihpb(nhpb)=ihpb(nhpb)+nres
1188 jhpb(nhpb)=jhpb(nhpb)+nres
1190 else if (restr_type.eq.10) then
1191 c Cross-lonk Markov-like potential
1192 call card_concat(controlcard)
1193 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1194 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1196 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1197 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1198 if (index(controlcard,"ZL").gt.0) then
1200 else if (index(controlcard,"ADH").gt.0) then
1202 else if (index(controlcard,"PDH").gt.0) then
1204 else if (index(controlcard,"DSS").gt.0) then
1209 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1210 & xlink(1,link_type))
1211 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1212 & xlink(2,link_type))
1213 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1214 & xlink(3,link_type))
1215 call reada(controlcard,"SIGMA",forcon(nhpb+1),
1216 & xlink(4,link_type))
1217 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1218 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1219 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1220 if (forcon(nhpb+1).le.0.0d0 .or.
1221 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1223 irestr_type(nhpb)=10
1224 c if (ibecarb(nhpb).gt.0) then
1225 c ihpb(nhpb)=ihpb(nhpb)+nres
1226 c jhpb(nhpb)=jhpb(nhpb)+nres
1228 if (ibecarb(nhpb).eq.3) then
1229 jhpb(nhpb)=jhpb(nhpb)+nres
1230 else if (ibecarb(nhpb).eq.2) then
1231 ihpb(nhpb)=ihpb(nhpb)+nres
1232 else if (ibecarb(nhpb).eq.1) then
1233 ihpb(nhpb)=ihpb(nhpb)+nres
1234 jhpb(nhpb)=jhpb(nhpb)+nres
1236 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1237 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1238 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1242 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1243 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1244 if (forcon(nhpb+1).gt.0.0d0) then
1246 if (dhpb1(nhpb).eq.0.0d0) then
1251 c if (ibecarb(nhpb).gt.0) then
1252 c ihpb(nhpb)=ihpb(nhpb)+nres
1253 c jhpb(nhpb)=jhpb(nhpb)+nres
1255 if (ibecarb(nhpb).eq.3) then
1256 jhpb(nhpb)=jhpb(nhpb)+nres
1257 else if (ibecarb(nhpb).eq.2) then
1258 ihpb(nhpb)=ihpb(nhpb)+nres
1259 else if (ibecarb(nhpb).eq.1) then
1260 ihpb(nhpb)=ihpb(nhpb)+nres
1261 jhpb(nhpb)=jhpb(nhpb)+nres
1263 if (dhpb(nhpb).eq.0.0d0)
1264 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1266 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1267 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1269 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1270 C if (forcon(nhpb+1).gt.0.0d0) then
1272 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1275 if (restr_type.eq.4) then
1276 write (iout,*) "The BFAC array"
1278 write (iout,'(i5,f10.5)') i,bfac(i)
1281 if (itype(i).eq.ntyp1) cycle
1283 if (itype(j).eq.ntyp1) cycle
1284 if (itype(i).eq.10) then
1289 if (itype(j).eq.10) then
1299 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1302 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1303 ihpb(nhpb)=i+nres*ii
1304 jhpb(nhpb)=j+nres*jj
1305 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1306 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1307 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1308 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1317 if (restr_type.eq.5) then
1318 restr_on_coord=.true.
1320 if (itype(i).eq.ntyp1) cycle
1321 bfac(i)=(scal_bfac/bfac(i))**2
1330 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1331 & fordepthmax=fordepth(i)
1334 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1337 if (nhpb.gt.nss) then
1338 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1339 & "The following",nhpb-nss,
1340 & " distance restraints have been imposed:",
1341 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1344 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1345 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1352 11 write (iout,*)"read_dist_restr: error reading reference structure"
1355 c-------------------------------------------------------------------------------
1356 subroutine read_saxs_constr
1357 implicit real*8 (a-h,o-z)
1358 include 'DIMENSIONS'
1362 include 'COMMON.CONTROL'
1363 include 'COMMON.CHAIN'
1364 include 'COMMON.IOUNITS'
1365 include 'COMMON.SBRIDGE'
1366 include 'COMMON.SAXS'
1367 double precision cm(3)
1369 write (iout,*) "Calling read_saxs nsaxs",nsaxs
1371 if (saxs_mode.eq.0) then
1372 c SAXS distance distribution
1374 read(inp,*) distsaxs(i),Psaxs(i)
1378 Cnorm = Cnorm + Psaxs(i)
1380 write (iout,*) "Cnorm",Cnorm
1382 Psaxs(i)=Psaxs(i)/Cnorm
1384 write (iout,*) "Normalized distance distribution from SAXS"
1386 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1390 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1392 write (iout,*) "Wsaxs0",Wsaxs0
1396 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1403 cm(j)=cm(j)+Csaxs(j,i)
1411 Csaxs(j,i)=Csaxs(j,i)-cm(j)
1414 write (iout,*) "SAXS sphere coordinates"
1416 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)