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,15.0d0)
43 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
44 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
45 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
46 if (lipthick.gt.0.0d0) then
47 bordliptop=(boxzsize+lipthick)/2.0
48 bordlipbot=bordliptop-lipthick
50 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
51 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
52 buflipbot=bordlipbot+lipbufthick
53 bufliptop=bordliptop-lipbufthick
54 if ((lipbufthick*2.0d0).gt.lipthick)
55 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
57 write(iout,*) "bordliptop=",bordliptop
58 write(iout,*) "bordlipbot=",bordlipbot
59 write(iout,*) "bufliptop=",bufliptop
60 write(iout,*) "buflipbot=",buflipbot
62 call readi(controlcard,'SHIELD',shield_mode,0)
63 write (iout,*) "SHIELD MODE",shield_mode
64 if (shield_mode.gt.0) then
65 pdbref=(index(controlcard,'PDBREF').gt.0)
66 if (index(controlcard,"CASC").gt.0) then
68 else if (index(controlcard,"SCONLY").gt.0) then
74 C VSolvSphere the volume of solving sphere
76 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
77 C there will be no distinction between proline peptide group and normal peptide
78 C group in case of shielding parameters
79 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
80 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
81 write (iout,*) VSolvSphere,VSolvSphere_div
82 C long axis of side chain
84 long_r_sidechain(i)=vbldsc0(1,i)
85 short_r_sidechain(i)=sigma0(i)
89 call readi(controlcard,'PDBOUT',outpdb,0)
90 call readi(controlcard,'MOL2OUT',outmol2,0)
91 refstr=(index(controlcard,'REFSTR').gt.0)
92 pdbref=(index(controlcard,'PDBREF').gt.0)
93 refstr = refstr .or. pdbref
94 write (iout,*) "REFSTR",refstr," PDBREF",pdbref
95 iscode=index(controlcard,'ONE_LETTER')
96 tree=(index(controlcard,'MAKE_TREE').gt.0)
97 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 write (iout,*) "with_dihed_constr ",with_dihed_constr,
100 & " CONSTR_DIST",constr_dist
101 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
102 write (iout,*) "with_theta_constr ",with_theta_constr
104 min_var=(index(controlcard,'MINVAR').gt.0)
105 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
106 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
107 print_fittest=(index(controlcard,'PRINT_FITTEST').gt.0)
108 call readi(controlcard,'NCUT',ncut,0)
110 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
113 call readi(controlcard,'NCLUST',nclust,5)
115 call readi(controlcard,'SYM',symetr,1)
116 write (iout,*) 'sym', symetr
117 call readi(controlcard,'NSTART',nstart,0)
118 call readi(controlcard,'NEND',nend,0)
119 call reada(controlcard,'ECUT',ecut,10.0d0)
120 call reada(controlcard,'PROB',prob_limit,0.99d0)
121 write (iout,*) "Probability limit",prob_limit
122 lgrp=(index(controlcard,'LGRP').gt.0)
123 caonly=(index(controlcard,'CA_ONLY').gt.0)
124 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
125 call readi(controlcard,'IOPT',iopt,2)
126 lefree = index(controlcard,"EFREE").gt.0
127 call readi(controlcard,'NTEMP',nT,1)
128 write (iout,*) "nT",nT
129 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
130 write (iout,*) "nT",nT
131 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
133 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
135 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
136 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
137 lprint_int=index(controlcard,"PRINT_INT") .gt.0
138 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
139 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
140 c if (constr_homology) tole=dmax1(tole,1.5d0)
141 write (iout,*) "with_homology_constr ",with_dihed_constr,
142 & " CONSTR_HOMOLOGY",constr_homology
143 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
144 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
145 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
146 write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD
147 call readi(controlcard,'NSAXS',nsaxs,0)
148 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
149 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
150 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
151 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
152 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
156 c--------------------------------------------------------------------------
159 C Read molecular data.
163 include 'COMMON.IOUNITS'
166 include 'COMMON.INTERACT'
167 include 'COMMON.LOCAL'
168 include 'COMMON.NAMES'
169 include 'COMMON.CHAIN'
170 include 'COMMON.FFIELD'
171 include 'COMMON.SBRIDGE'
172 include 'COMMON.HEADER'
173 include 'COMMON.CONTROL'
174 include 'COMMON.CONTACTS'
175 include 'COMMON.TIME1'
176 include 'COMMON.TORCNSTR'
177 include 'COMMON.SHIELD'
178 include 'COMMON.SAXS'
180 include 'COMMON.INFO'
182 character*4 sequence(maxres)
183 character*800 weightcard,controlcard
185 double precision x(maxvar)
186 double precision phihel,phibet,sigmahel,sigmabet,sumv,
188 integer itype_pdb(maxres)
190 integer i,j,kkk,i1,i2,it1,it2,tperm,ii,iperm
194 C Read weights of the subsequent energy terms.
195 call card_concat(weightcard)
196 call reada(weightcard,'WSC',wsc,1.0d0)
197 call reada(weightcard,'WLONG',wsc,wsc)
198 call reada(weightcard,'WSCP',wscp,1.0d0)
199 call reada(weightcard,'WELEC',welec,1.0D0)
200 call reada(weightcard,'WVDWPP',wvdwpp,welec)
201 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
202 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
203 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
204 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
205 call reada(weightcard,'WTURN3',wturn3,1.0D0)
206 call reada(weightcard,'WTURN4',wturn4,1.0D0)
207 call reada(weightcard,'WTURN6',wturn6,1.0D0)
208 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
209 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
210 call reada(weightcard,'WBOND',wbond,1.0D0)
211 call reada(weightcard,'WTOR',wtor,1.0D0)
212 call reada(weightcard,'WTORD',wtor_d,1.0D0)
213 call reada(weightcard,'WANG',wang,1.0D0)
214 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
215 call reada(weightcard,'WSAXS',wsaxs,0.0D0)
216 call reada(weightcard,'SCAL14',scal14,0.4D0)
217 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
218 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
219 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
220 if (index(weightcard,'SOFT').gt.0) ipot=6
221 call reada(weightcard,"D0CM",d0cm,3.78d0)
222 call reada(weightcard,"AKCM",akcm,15.1d0)
223 call reada(weightcard,"AKTH",akth,11.0d0)
224 call reada(weightcard,"AKCT",akct,12.0d0)
225 call reada(weightcard,"V1SS",v1ss,-1.08d0)
226 call reada(weightcard,"V2SS",v2ss,7.61d0)
227 call reada(weightcard,"V3SS",v3ss,13.7d0)
228 call reada(weightcard,"EBR",ebr,-5.50D0)
229 call reada(weightcard,'WSHIELD',wshield,1.0d0)
230 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
231 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
232 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
233 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
234 call reada(weightcard,'WLT',wliptran,0.0D0)
235 call reada(weightcard,"ATRISS",atriss,0.301D0)
236 call reada(weightcard,"BTRISS",btriss,0.021D0)
237 call reada(weightcard,"CTRISS",ctriss,1.001D0)
238 call reada(weightcard,"DTRISS",dtriss,1.001D0)
239 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
241 dyn_ss_mask(i)=.false.
245 dyn_ssbond_ij(i,j)=1.0d300
248 call reada(weightcard,"HT",Ht,0.0D0)
250 ss_depth=ebr/wsc-0.25*eps(1,1)
251 Ht=Ht/wsc-0.25*eps(1,1)
252 akcm=akcm*wstrain/wsc
253 akth=akth*wstrain/wsc
254 akct=akct*wstrain/wsc
255 v1ss=v1ss*wstrain/wsc
256 v2ss=v2ss*wstrain/wsc
257 v3ss=v3ss*wstrain/wsc
259 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
261 write (iout,'(/a)') "Disulfide bridge parameters:"
262 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
263 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
264 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
265 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
267 write (iout,*) "Parameters of the 'trisulfide' potential"
268 write (iout,*) "ATRISS=", atriss
269 write (iout,*) "BTRISS=", btriss
270 write (iout,*) "CTRISS=", ctriss
271 write (iout,*) "DTRISS=", dtriss
273 C 12/1/95 Added weight for the multi-body term WCORR
274 call reada(weightcard,'WCORRH',wcorr,1.0D0)
275 if (wcorr4.gt.0.0d0) wcorr=wcorr4
295 weights(28)=wdfa_dist
298 weights(31)=wdfa_beta
299 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
300 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
301 & wturn4,wturn6,wsccor
302 10 format (/'Energy-term weights (unscaled):'//
303 & 'WSCC= ',f10.6,' (SC-SC)'/
304 & 'WSCP= ',f10.6,' (SC-p)'/
305 & 'WELEC= ',f10.6,' (p-p electr)'/
306 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
307 & 'WBOND= ',f10.6,' (stretching)'/
308 & 'WANG= ',f10.6,' (bending)'/
309 & 'WSCLOC= ',f10.6,' (SC local)'/
310 & 'WTOR= ',f10.6,' (torsional)'/
311 & 'WTORD= ',f10.6,' (double torsional)'/
312 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
313 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
314 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
315 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
316 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
317 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
318 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
319 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
320 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
322 if (wcorr4.gt.0.0d0) then
323 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
324 & 'between contact pairs of peptide groups'
325 write (iout,'(2(a,f5.3/))')
326 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
327 & 'Range of quenching the correlation terms:',2*delt_corr
328 else if (wcorr.gt.0.0d0) then
329 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
330 & 'between contact pairs of peptide groups'
332 write (iout,'(a,f8.3)')
333 & 'Scaling factor of 1,4 SC-p interactions:',scal14
334 write (iout,'(a,f8.3)')
335 & 'General scaling factor of SC-p interactions:',scalscp
336 r0_corr=cutoff_corr-delt_corr
338 aad(i,1)=scalscp*aad(i,1)
339 aad(i,2)=scalscp*aad(i,2)
340 bad(i,1)=scalscp*bad(i,1)
341 bad(i,2)=scalscp*bad(i,2)
344 write (iout,'(/a/)') "DFA pseudopotential parameters:"
345 write (iout,'(a,f10.6,a)')
346 & "WDFAD= ",wdfa_dist," (distance)",
347 & "WDFAT= ",wdfa_tor," (backbone angles)",
348 & "WDFAN= ",wdfa_nei," (neighbors)",
349 & "WDFAB= ",wdfa_beta," (beta structure)"
352 c print *,'indpdb=',indpdb,' pdbref=',pdbref
354 C Read sequence if not taken from the pdb file.
355 if (iscode.gt.0) then
356 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
358 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
360 C Convert sequence to numeric code
362 itype(i)=rescode(i,sequence(i),iscode)
365 c print '(20i4)',(itype(i),i=1,nres)
369 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
371 if (itype(i).eq.ntyp1) then
375 else if (iabs(itype(i+1)).ne.20) then
377 else if (iabs(itype(i)).ne.20) then
384 write (iout,*) "ITEL"
386 write (iout,*) i,itype(i),itel(i)
389 c print *,'Call Read_Bridge.'
391 C this fragment reads diheadral constrains
394 c print *,'NNT=',NNT,' NCT=',NCT
395 call seq2chains(nres,itype,nchain,chain_length,chain_border,
397 write(iout,*) "nres",nres," nchain",nchain
399 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
402 call chain_symmetry(nchain,nres,itype,chain_border,
403 & chain_length,npermchain,tabpermchain)
405 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
408 write(iout,*) "residue permutations"
410 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
412 if (itype(1).eq.ntyp1) nnt=2
413 if (itype(nres).eq.ntyp1) nct=nct-1
414 if (nstart.lt.nnt) nstart=nnt
415 if (nend.gt.nct .or. nend.eq.0) nend=nct
416 write (iout,*) "nstart",nstart," nend",nend
419 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
420 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
422 print*, 'init_dfa_vars finished!'
424 print*, 'read_dfa_info finished!'
427 if (with_dihed_constr) then
429 read (inp,*) ndih_constr
430 if (ndih_constr.gt.0) then
433 C write (iout,*) 'FTORS',ftors
434 C ftors is the force constant for torsional quartic constrains
435 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
438 & 'There are',ndih_constr,' constraints on phi angles.'
440 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
444 phi0(i)=deg2rad*phi0(i)
445 drange(i)=deg2rad*drange(i)
447 else if (ndih_constr.lt.0) then
449 call card_concat(controlcard)
450 call reada(controlcard,"PHIHEL",phihel,50.0D0)
451 call reada(controlcard,"PHIBET",phibet,180.0D0)
452 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
453 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
454 call reada(controlcard,"WDIHC",wdihc,0.591d0)
455 write (iout,*) "Weight of the dihedral restraint term",wdihc
456 read(inp,'(9x,3f7.3)')
457 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
458 write (iout,*) "The secprob array"
460 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
464 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
465 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
466 ndih_constr=ndih_constr+1
467 idih_constr(ndih_constr)=i
470 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
471 sumv=sumv+vpsipred(j,ndih_constr)
474 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
476 phibound(1,ndih_constr)=phihel*deg2rad
477 phibound(2,ndih_constr)=phibet*deg2rad
478 sdihed(1,ndih_constr)=sigmahel*deg2rad
479 sdihed(2,ndih_constr)=sigmabet*deg2rad
483 & 'There are',ndih_constr,
484 & ' bimodal restraints on gamma angles.'
486 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
487 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
488 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
489 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
490 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
491 & (vpsipred(j,i),j=1,3)
494 endif ! endif ndif_constr.gt.0
495 endif ! with_dihed_constr
496 if (with_theta_constr) then
497 C with_theta_constr is keyword allowing for occurance of theta constrains
498 read (inp,*) ntheta_constr
499 C ntheta_constr is the number of theta constrains
500 if (ntheta_constr.gt.0) then
502 read (inp,*) (itheta_constr(i),theta_constr0(i),
503 & theta_drange(i),for_thet_constr(i),
505 C the above code reads from 1 to ntheta_constr
506 C itheta_constr(i) residue i for which is theta_constr
507 C theta_constr0 the global minimum value
508 C theta_drange is range for which there is no energy penalty
509 C for_thet_constr is the force constant for quartic energy penalty
512 & 'There are',ntheta_constr,' constraints on phi angles.'
514 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
520 theta_constr0(i)=deg2rad*theta_constr0(i)
521 theta_drange(i)=deg2rad*theta_drange(i)
523 C if(me.eq.king.or..not.out1file)
524 C & write (iout,*) 'FTORS',ftors
525 C do i=1,ntheta_constr
526 C ii = itheta_constr(i)
527 C thetabound(1,ii) = phi0(i)-drange(i)
528 C thetabound(2,ii) = phi0(i)+drange(i)
530 endif ! ntheta_constr.gt.0
531 endif! with_theta_constr
532 if (constr_homology.gt.0) then
533 c write (iout,*) "About to call read_constr_homology"
535 call read_constr_homology
536 c write (iout,*) "Exit read_constr_homology"
538 if (indpdb.gt.0 .or. pdbref) then
542 cref(j,i)=crefjlee(j,i)
547 write (iout,*) "Array C"
549 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
550 & (c(j,i+nres),j=1,3)
552 write (iout,*) "Array Cref"
554 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
555 & (cref(j,i+nres),j=1,3)
559 call int_from_cart1(.false.)
560 call sc_loc_geom(.false.)
564 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
568 dc(j,i)=c(j,i+1)-c(j,i)
569 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
574 dc(j,i+nres)=c(j,i+nres)-c(j,i)
575 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
582 write (iout,*) "calling read_saxs_consrtr",nsaxs
583 if (nsaxs.gt.0) call read_saxs_constr
586 c read(inp,'(a)') pdbfile
587 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
588 c open(ipdbin,file=pdbfile,status='old',err=33)
590 c 33 write (iout,'(a)') 'Error opening PDB file.'
593 c print *,'Begin reading pdb data'
595 c print *,'Finished reading pdb data'
596 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
598 c itype_pdb(i)=itype(i)
601 c write (iout,'(a,i3)') 'nsup=',nsup
603 c if (nsup.le.(nct-nnt+1)) then
604 c do i=0,nct-nnt+1-nsup
605 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
611 c & 'Error - sequences to be superposed do not match.'
614 c do i=0,nsup-(nct-nnt+1)
615 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
617 c nstart_sup=nstart_sup+i
623 c & 'Error - sequences to be superposed do not match.'
626 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
627 c & ' nstart_seq=',nstart_seq
631 write (iout,*) "molread: REFSTR",refstr
633 if (.not.pdbref) then
634 call read_angles(inp,*38)
636 38 write (iout,'(a)') 'Error reading reference structure.'
638 call mp_stopall(Error_Msg)
640 stop 'Error reading reference structure'
652 c call contact(.true.,ncont_ref,icont_ref)
655 C write (iout,'(/a,i3,a)')
656 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
657 write (iout,'(20i4)') (iss(i),i=1,ns)
659 write(iout,*)"Running with dynamic disulfide-bond formation"
661 write (iout,'(/a/)') 'Pre-formed links are:'
667 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
668 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
674 if (ns.gt.0.and.dyn_ss) then
678 forcon(i-nss)=forcon(i)
685 dyn_ss_mask(iss(i))=.true.
688 c Read distance restraints
689 if (constr_dist.gt.0) then
690 call read_dist_constr
695 c-----------------------------------------------------------------------------
696 logical function seq_comp(itypea,itypeb,length)
698 integer length,itypea(length),itypeb(length)
701 if (itypea(i).ne.itypeb(i)) then
709 c-----------------------------------------------------------------------------
710 subroutine read_bridge
711 C Read information about disulfide bridges.
714 include 'COMMON.IOUNITS'
717 include 'COMMON.INTERACT'
718 include 'COMMON.LOCAL'
719 include 'COMMON.NAMES'
720 include 'COMMON.CHAIN'
721 include 'COMMON.FFIELD'
722 include 'COMMON.SBRIDGE'
723 include 'COMMON.HEADER'
724 include 'COMMON.CONTROL'
725 include 'COMMON.TIME1'
727 include 'COMMON.INFO'
730 C Read bridging residues.
731 read (inp,*) ns,(iss(i),i=1,ns)
733 C Check whether the specified bridging residues are cystines.
735 if (itype(iss(i)).ne.1) then
736 write (iout,'(2a,i3,a)')
737 & 'Do you REALLY think that the residue ',
738 & restyp(itype(iss(i))),i,
739 & ' can form a disulfide bridge?!!!'
740 write (*,'(2a,i3,a)')
741 & 'Do you REALLY think that the residue ',
742 & restyp(itype(iss(i))),i,
743 & ' can form a disulfide bridge?!!!'
745 call mp_stopall(error_msg)
751 C Read preformed bridges.
753 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
756 C Check if the residues involved in bridges are in the specified list of
760 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
761 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
762 write (iout,'(a,i3,a)') 'Disulfide pair',i,
763 & ' contains residues present in other pairs.'
764 write (*,'(a,i3,a)') 'Disulfide pair',i,
765 & ' contains residues present in other pairs.'
767 call mp_stopall(error_msg)
774 if (ihpb(i).eq.iss(j)) goto 10
776 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
779 if (jhpb(i).eq.iss(j)) goto 20
781 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
794 c----------------------------------------------------------------------------
795 subroutine read_angles(kanal,*)
800 include 'COMMON.CHAIN'
801 include 'COMMON.IOUNITS'
803 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
804 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
805 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
806 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
808 theta(i)=deg2rad*theta(i)
809 phi(i)=deg2rad*phi(i)
810 alph(i)=deg2rad*alph(i)
811 omeg(i)=deg2rad*omeg(i)
816 c----------------------------------------------------------------------------
817 subroutine reada(rekord,lancuch,wartosc,default)
819 character*(*) rekord,lancuch
820 double precision wartosc,default
823 iread=index(rekord,lancuch)
828 iread=iread+ilen(lancuch)+1
829 read (rekord(iread:),*) wartosc
832 c----------------------------------------------------------------------------
833 subroutine multreada(rekord,lancuch,tablica,dim,default)
836 double precision tablica(dim),default
837 character*(*) rekord,lancuch
843 iread=index(rekord,lancuch)
844 if (iread.eq.0) return
845 iread=iread+ilen(lancuch)+1
846 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
849 c----------------------------------------------------------------------------
850 subroutine readi(rekord,lancuch,wartosc,default)
852 character*(*) rekord,lancuch
853 integer wartosc,default
856 iread=index(rekord,lancuch)
861 iread=iread+ilen(lancuch)+1
862 read (rekord(iread:),*) wartosc
865 C----------------------------------------------------------------------
866 subroutine multreadi(rekord,lancuch,tablica,dim,default)
869 integer tablica(dim),default
870 character*(*) rekord,lancuch
877 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
878 if (iread.eq.0) return
879 iread=iread+ilen(lancuch)+1
880 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
884 c----------------------------------------------------------------------------
885 subroutine card_concat(card)
887 include 'COMMON.IOUNITS'
889 character*80 karta,ucase
891 read (inp,'(a)') karta
894 do while (karta(80:80).eq.'&')
895 card=card(:ilen(card)+1)//karta(:79)
896 read (inp,'(a)') karta
899 card=card(:ilen(card)+1)//karta
902 c----------------------------------------------------------------------------
911 include 'COMMON.IOUNITS'
912 include 'COMMON.CONTROL'
913 integer lenpre,lenpot,ilen
915 character*16 cformat,cprint
917 integer lenint,lenout
918 call getenv('INPUT',prefix)
919 call getenv('OUTPUT',prefout)
920 call getenv('INTIN',prefintin)
921 call getenv('COORD',cformat)
922 call getenv('PRINTCOOR',cprint)
923 call getenv('SCRATCHDIR',scratchdir)
926 if (index(ucase(cformat),'CX').gt.0) then
933 lenint=ilen(prefintin)
934 C Get the names and open the input files
935 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
937 write (liczba,'(bz,i3.3)') me
938 outname=prefout(:lenout)//'_clust.out_'//liczba
940 outname=prefout(:lenout)//'_clust.out'
943 intinname=prefintin(:lenint)//'.bx'
944 else if (from_cx) then
945 intinname=prefintin(:lenint)//'.cx'
947 intinname=prefintin(:lenint)//'.int'
949 rmsname=prefintin(:lenint)//'.rms'
950 open (jplot,file=prefout(:ilen(prefout))//'.tex',
952 open (jrms,file=rmsname,status='unknown')
953 open(iout,file=outname,status='unknown')
954 C Get parameter filenames and open the parameter files.
955 call getenv('BONDPAR',bondname)
956 open (ibond,file=bondname,status='old')
957 call getenv('THETPAR',thetname)
958 open (ithep,file=thetname,status='old')
959 call getenv('ROTPAR',rotname)
960 open (irotam,file=rotname,status='old')
961 call getenv('TORPAR',torname)
962 open (itorp,file=torname,status='old')
964 call getenv('TORDPAR',tordname)
965 open (itordp,file=tordname,status='old')
967 call getenv('FOURIER',fouriername)
968 open (ifourier,file=fouriername,status='old')
969 call getenv('ELEPAR',elename)
970 open (ielep,file=elename,status='old')
971 call getenv('SIDEPAR',sidename)
972 open (isidep,file=sidename,status='old')
973 call getenv('SIDEP',sidepname)
974 open (isidep1,file=sidepname,status="old")
975 call getenv('SCCORPAR',sccorname)
976 open (isccor,file=sccorname,status="old")
977 call getenv('LIPTRANPAR',liptranname)
978 open (iliptranpar,file=liptranname,status='old')
981 C 8/9/01 In the newest version SCp interaction constants are read from a file
982 C Use -DOLDSCP to use hard-coded constants instead.
984 call getenv('SCPPAR',scpname)
985 open (iscpp,file=scpname,status='old')
989 c--------------------------------------------------------------------------
990 subroutine read_dist_constr
991 implicit real*8 (a-h,o-z)
996 include 'COMMON.CONTROL'
997 include 'COMMON.CHAIN'
998 include 'COMMON.IOUNITS'
999 include 'COMMON.SBRIDGE'
1000 include 'COMMON.INTERACT'
1001 integer ifrag_(2,100),ipair_(2,100)
1002 double precision wfrag_(100),wpair_(100)
1003 character*500 controlcard
1004 logical lprn /.true./
1005 logical normalize,next
1007 double precision scal_bfac
1008 double precision xlink(4,0:4) /
1010 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
1011 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
1012 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
1013 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
1014 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
1015 write (iout,*) "Calling read_dist_constr"
1016 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1022 call card_concat(controlcard)
1023 next = index(controlcard,"NEXT").gt.0
1024 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
1025 write (iout,*) "restr_type",restr_type
1026 call readi(controlcard,"NFRAG",nfrag_,0)
1027 call readi(controlcard,"NFRAG",nfrag_,0)
1028 call readi(controlcard,"NPAIR",npair_,0)
1029 call readi(controlcard,"NDIST",ndist_,0)
1030 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1031 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
1032 if (restr_type.eq.10)
1033 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
1034 if (restr_type.eq.12)
1035 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
1036 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1037 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1038 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1039 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1040 normalize = index(controlcard,"NORMALIZE").gt.0
1041 write (iout,*) "WBOLTZD",wboltzd
1042 write (iout,*) "SCAL_PEAK",scal_peak
1043 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1044 write (iout,*) "IFRAG"
1046 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1048 write (iout,*) "IPAIR"
1050 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1052 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
1054 read(inp,'(a)') pdbfile
1056 & "Distance restraints will be constructed from structure ",pdbfile
1057 open(ipdbin,file=pdbfile,status='old',err=11)
1058 call readpdb(.true.)
1063 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1064 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1065 & ifrag_(2,i)=nstart_sup+nsup-1
1066 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1068 if (wfrag_(i).eq.0.0d0) cycle
1069 do j=ifrag_(1,i),ifrag_(2,i)-1
1070 do k=j+1,ifrag_(2,i)
1071 c write (iout,*) "j",j," k",k
1073 if (restr_type.eq.1) then
1079 forcon(nhpb)=wfrag_(i)
1080 else if (constr_dist.eq.2) then
1081 if (ddjk.le.dist_cut) then
1087 forcon(nhpb)=wfrag_(i)
1089 else if (restr_type.eq.3) then
1095 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1097 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
1098 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1103 if (wpair_(i).eq.0.0d0) cycle
1111 do j=ifrag_(1,ii),ifrag_(2,ii)
1112 do k=ifrag_(1,jj),ifrag_(2,jj)
1114 if (restr_type.eq.1) then
1120 forcon(nhpb)=wpair_(i)
1121 else if (constr_dist.eq.2) then
1122 if (ddjk.le.dist_cut) then
1128 forcon(nhpb)=wpair_(i)
1130 else if (restr_type.eq.3) then
1136 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1138 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1139 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1145 write (iout,*) "Distance restraints as read from input"
1147 if (restr_type.eq.12) then
1148 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1149 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1150 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1151 & fordepth_peak(nhpb_peak+1),npeak
1152 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1153 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1154 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1155 c & fordepth_peak(nhpb_peak+1),npeak
1156 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1157 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1158 nhpb_peak=nhpb_peak+1
1159 irestr_type_peak(nhpb_peak)=12
1160 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1162 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1163 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1164 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1165 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1166 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1167 if (ibecarb_peak(nhpb_peak).eq.3) then
1168 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1169 else if (ibecarb_peak(nhpb_peak).eq.2) then
1170 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1171 else if (ibecarb_peak(nhpb_peak).eq.1) then
1172 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1173 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1175 else if (restr_type.eq.11) then
1176 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1177 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1178 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1179 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1181 irestr_type(nhpb)=11
1182 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1183 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1184 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1185 c if (ibecarb(nhpb).gt.0) then
1186 c ihpb(nhpb)=ihpb(nhpb)+nres
1187 c jhpb(nhpb)=jhpb(nhpb)+nres
1189 if (ibecarb(nhpb).eq.3) then
1190 ihpb(nhpb)=ihpb(nhpb)+nres
1191 else if (ibecarb(nhpb).eq.2) then
1192 ihpb(nhpb)=ihpb(nhpb)+nres
1193 else if (ibecarb(nhpb).eq.1) then
1194 ihpb(nhpb)=ihpb(nhpb)+nres
1195 jhpb(nhpb)=jhpb(nhpb)+nres
1197 else if (restr_type.eq.10) then
1198 c Cross-lonk Markov-like potential
1199 call card_concat(controlcard)
1200 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1201 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1203 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1204 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1205 if (index(controlcard,"ZL").gt.0) then
1207 else if (index(controlcard,"ADH").gt.0) then
1209 else if (index(controlcard,"PDH").gt.0) then
1211 else if (index(controlcard,"DSS").gt.0) then
1216 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1217 & xlink(1,link_type))
1218 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1219 & xlink(2,link_type))
1220 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1221 & xlink(3,link_type))
1222 call reada(controlcard,"SIGMA",forcon(nhpb+1),
1223 & xlink(4,link_type))
1224 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1225 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1226 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1227 if (forcon(nhpb+1).le.0.0d0 .or.
1228 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1230 irestr_type(nhpb)=10
1231 c if (ibecarb(nhpb).gt.0) then
1232 c ihpb(nhpb)=ihpb(nhpb)+nres
1233 c jhpb(nhpb)=jhpb(nhpb)+nres
1235 if (ibecarb(nhpb).eq.3) then
1236 jhpb(nhpb)=jhpb(nhpb)+nres
1237 else if (ibecarb(nhpb).eq.2) then
1238 ihpb(nhpb)=ihpb(nhpb)+nres
1239 else if (ibecarb(nhpb).eq.1) then
1240 ihpb(nhpb)=ihpb(nhpb)+nres
1241 jhpb(nhpb)=jhpb(nhpb)+nres
1243 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1244 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1245 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1249 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1250 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1251 if (forcon(nhpb+1).gt.0.0d0) then
1253 if (dhpb1(nhpb).eq.0.0d0) then
1258 c if (ibecarb(nhpb).gt.0) then
1259 c ihpb(nhpb)=ihpb(nhpb)+nres
1260 c jhpb(nhpb)=jhpb(nhpb)+nres
1262 if (ibecarb(nhpb).eq.3) then
1263 jhpb(nhpb)=jhpb(nhpb)+nres
1264 else if (ibecarb(nhpb).eq.2) then
1265 ihpb(nhpb)=ihpb(nhpb)+nres
1266 else if (ibecarb(nhpb).eq.1) then
1267 ihpb(nhpb)=ihpb(nhpb)+nres
1268 jhpb(nhpb)=jhpb(nhpb)+nres
1270 if (dhpb(nhpb).eq.0.0d0)
1271 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1273 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1274 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1276 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1277 C if (forcon(nhpb+1).gt.0.0d0) then
1279 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1282 if (restr_type.eq.4) then
1283 write (iout,*) "The BFAC array"
1285 write (iout,'(i5,f10.5)') i,bfac(i)
1288 if (itype(i).eq.ntyp1) cycle
1290 if (itype(j).eq.ntyp1) cycle
1291 if (itype(i).eq.10) then
1296 if (itype(j).eq.10) then
1306 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1309 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1310 ihpb(nhpb)=i+nres*ii
1311 jhpb(nhpb)=j+nres*jj
1312 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1313 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1314 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1315 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1324 if (restr_type.eq.5) then
1325 restr_on_coord=.true.
1327 if (itype(i).eq.ntyp1) cycle
1328 bfac(i)=(scal_bfac/bfac(i))**2
1337 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1338 & fordepthmax=fordepth(i)
1341 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1344 if (nhpb.gt.nss) then
1345 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1346 & "The following",nhpb-nss,
1347 & " distance restraints have been imposed:",
1348 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1351 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1352 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1359 11 write (iout,*)"read_dist_restr: error reading reference structure"
1362 c-------------------------------------------------------------------------------
1363 subroutine read_saxs_constr
1364 implicit real*8 (a-h,o-z)
1365 include 'DIMENSIONS'
1369 include 'COMMON.CONTROL'
1370 include 'COMMON.CHAIN'
1371 include 'COMMON.IOUNITS'
1372 include 'COMMON.SBRIDGE'
1373 include 'COMMON.SAXS'
1374 double precision cm(3)
1376 write (iout,*) "Calling read_saxs nsaxs",nsaxs
1378 if (saxs_mode.eq.0) then
1379 c SAXS distance distribution
1381 read(inp,*) distsaxs(i),Psaxs(i)
1385 Cnorm = Cnorm + Psaxs(i)
1387 write (iout,*) "Cnorm",Cnorm
1389 Psaxs(i)=Psaxs(i)/Cnorm
1391 write (iout,*) "Normalized distance distribution from SAXS"
1393 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1397 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1399 write (iout,*) "Wsaxs0",Wsaxs0
1403 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1410 cm(j)=cm(j)+Csaxs(j,i)
1418 Csaxs(j,i)=Csaxs(j,i)-cm(j)
1421 write (iout,*) "SAXS sphere coordinates"
1423 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)