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,"ATRISS",atriss,0.301D0)
238 call reada(weightcard,"BTRISS",btriss,0.021D0)
239 call reada(weightcard,"CTRISS",ctriss,1.001D0)
240 call reada(weightcard,"DTRISS",dtriss,1.001D0)
241 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
243 dyn_ss_mask(i)=.false.
247 dyn_ssbond_ij(i,j)=1.0d300
250 call reada(weightcard,"HT",Ht,0.0D0)
252 ss_depth=ebr/wsc-0.25*eps(1,1)
253 Ht=Ht/wsc-0.25*eps(1,1)
254 akcm=akcm*wstrain/wsc
255 akth=akth*wstrain/wsc
256 akct=akct*wstrain/wsc
257 v1ss=v1ss*wstrain/wsc
258 v2ss=v2ss*wstrain/wsc
259 v3ss=v3ss*wstrain/wsc
261 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
263 write (iout,'(/a)') "Disulfide bridge parameters:"
264 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
265 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
266 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
267 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
269 write (iout,*) "Parameters of the 'trisulfide' potential"
270 write (iout,*) "ATRISS=", atriss
271 write (iout,*) "BTRISS=", btriss
272 write (iout,*) "CTRISS=", ctriss
273 write (iout,*) "DTRISS=", dtriss
275 C 12/1/95 Added weight for the multi-body term WCORR
276 call reada(weightcard,'WCORRH',wcorr,1.0D0)
277 if (wcorr4.gt.0.0d0) wcorr=wcorr4
297 weights(28)=wdfa_dist
300 weights(31)=wdfa_beta
301 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
302 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
303 & wturn4,wturn6,wsccor
304 10 format (/'Energy-term weights (unscaled):'//
305 & 'WSCC= ',f10.6,' (SC-SC)'/
306 & 'WSCP= ',f10.6,' (SC-p)'/
307 & 'WELEC= ',f10.6,' (p-p electr)'/
308 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
309 & 'WBOND= ',f10.6,' (stretching)'/
310 & 'WANG= ',f10.6,' (bending)'/
311 & 'WSCLOC= ',f10.6,' (SC local)'/
312 & 'WTOR= ',f10.6,' (torsional)'/
313 & 'WTORD= ',f10.6,' (double torsional)'/
314 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
315 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
316 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
317 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
318 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
319 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
320 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
321 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
322 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
324 if (wcorr4.gt.0.0d0) then
325 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
326 & 'between contact pairs of peptide groups'
327 write (iout,'(2(a,f5.3/))')
328 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
329 & 'Range of quenching the correlation terms:',2*delt_corr
330 else if (wcorr.gt.0.0d0) then
331 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
332 & 'between contact pairs of peptide groups'
334 write (iout,'(a,f8.3)')
335 & 'Scaling factor of 1,4 SC-p interactions:',scal14
336 write (iout,'(a,f8.3)')
337 & 'General scaling factor of SC-p interactions:',scalscp
338 r0_corr=cutoff_corr-delt_corr
340 aad(i,1)=scalscp*aad(i,1)
341 aad(i,2)=scalscp*aad(i,2)
342 bad(i,1)=scalscp*bad(i,1)
343 bad(i,2)=scalscp*bad(i,2)
346 write (iout,'(/a/)') "DFA pseudopotential parameters:"
347 write (iout,'(a,f10.6,a)')
348 & "WDFAD= ",wdfa_dist," (distance)",
349 & "WDFAT= ",wdfa_tor," (backbone angles)",
350 & "WDFAN= ",wdfa_nei," (neighbors)",
351 & "WDFAB= ",wdfa_beta," (beta structure)"
354 c print *,'indpdb=',indpdb,' pdbref=',pdbref
356 C Read sequence if not taken from the pdb file.
357 if (iscode.gt.0) then
358 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
360 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
362 C Convert sequence to numeric code
364 itype(i)=rescode(i,sequence(i),iscode)
367 c print '(20i4)',(itype(i),i=1,nres)
371 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
373 if (itype(i).eq.ntyp1) then
377 else if (iabs(itype(i+1)).ne.20) then
379 else if (iabs(itype(i)).ne.20) then
386 write (iout,*) "ITEL"
388 write (iout,*) i,itype(i),itel(i)
391 c print *,'Call Read_Bridge.'
393 C this fragment reads diheadral constrains
396 c print *,'NNT=',NNT,' NCT=',NCT
397 call seq2chains(nres,itype,nchain,chain_length,chain_border,
399 write(iout,*) "nres",nres," nchain",nchain
401 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
404 call chain_symmetry(nchain,nres,itype,chain_border,
405 & chain_length,npermchain,tabpermchain)
407 write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
410 write(iout,*) "residue permutations"
412 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
414 if (itype(1).eq.ntyp1) nnt=2
415 if (itype(nres).eq.ntyp1) nct=nct-1
416 if (nstart.lt.nnt) nstart=nnt
417 if (nend.gt.nct .or. nend.eq.0) nend=nct
418 write (iout,*) "nstart",nstart," nend",nend
421 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
422 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
424 print*, 'init_dfa_vars finished!'
426 print*, 'read_dfa_info finished!'
429 C If the reference structure is not read set the superposition
436 if (with_dihed_constr) then
438 read (inp,*) ndih_constr
439 if (ndih_constr.gt.0) then
442 C write (iout,*) 'FTORS',ftors
443 C ftors is the force constant for torsional quartic constrains
444 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
447 & 'There are',ndih_constr,' constraints on phi angles.'
449 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
453 phi0(i)=deg2rad*phi0(i)
454 drange(i)=deg2rad*drange(i)
456 else if (ndih_constr.lt.0) then
458 call card_concat(controlcard)
459 call reada(controlcard,"PHIHEL",phihel,50.0D0)
460 call reada(controlcard,"PHIBET",phibet,180.0D0)
461 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
462 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
463 call reada(controlcard,"WDIHC",wdihc,0.591d0)
464 write (iout,*) "Weight of the dihedral restraint term",wdihc
465 read(inp,'(9x,3f7.3)')
466 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
467 write (iout,*) "The secprob array"
469 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
473 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
474 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
475 ndih_constr=ndih_constr+1
476 idih_constr(ndih_constr)=i
479 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
480 sumv=sumv+vpsipred(j,ndih_constr)
483 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
485 phibound(1,ndih_constr)=phihel*deg2rad
486 phibound(2,ndih_constr)=phibet*deg2rad
487 sdihed(1,ndih_constr)=sigmahel*deg2rad
488 sdihed(2,ndih_constr)=sigmabet*deg2rad
492 & 'There are',ndih_constr,
493 & ' bimodal restraints on gamma angles.'
495 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
496 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
497 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
498 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
499 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
500 & (vpsipred(j,i),j=1,3)
503 endif ! endif ndif_constr.gt.0
504 endif ! with_dihed_constr
505 if (with_theta_constr) then
506 C with_theta_constr is keyword allowing for occurance of theta constrains
507 read (inp,*) ntheta_constr
508 C ntheta_constr is the number of theta constrains
509 if (ntheta_constr.gt.0) then
511 read (inp,*) (itheta_constr(i),theta_constr0(i),
512 & theta_drange(i),for_thet_constr(i),
514 C the above code reads from 1 to ntheta_constr
515 C itheta_constr(i) residue i for which is theta_constr
516 C theta_constr0 the global minimum value
517 C theta_drange is range for which there is no energy penalty
518 C for_thet_constr is the force constant for quartic energy penalty
521 & 'There are',ntheta_constr,' constraints on phi angles.'
523 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
529 theta_constr0(i)=deg2rad*theta_constr0(i)
530 theta_drange(i)=deg2rad*theta_drange(i)
532 C if(me.eq.king.or..not.out1file)
533 C & write (iout,*) 'FTORS',ftors
534 C do i=1,ntheta_constr
535 C ii = itheta_constr(i)
536 C thetabound(1,ii) = phi0(i)-drange(i)
537 C thetabound(2,ii) = phi0(i)+drange(i)
539 endif ! ntheta_constr.gt.0
540 endif! with_theta_constr
541 if (constr_homology.gt.0) then
542 c write (iout,*) "About to call read_constr_homology"
544 call read_constr_homology
545 c write (iout,*) "Exit read_constr_homology"
547 if (indpdb.gt.0 .or. pdbref) then
551 cref(j,i)=crefjlee(j,i)
556 write (iout,*) "Array C"
558 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
559 & (c(j,i+nres),j=1,3)
561 write (iout,*) "Array Cref"
563 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
564 & (cref(j,i+nres),j=1,3)
568 call int_from_cart1(.false.)
569 call sc_loc_geom(.false.)
573 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
577 dc(j,i)=c(j,i+1)-c(j,i)
578 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
583 dc(j,i+nres)=c(j,i+nres)-c(j,i)
584 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
591 write (iout,*) "calling read_saxs_consrtr",nsaxs
592 if (nsaxs.gt.0) call read_saxs_constr
595 c read(inp,'(a)') pdbfile
596 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
597 c open(ipdbin,file=pdbfile,status='old',err=33)
599 c 33 write (iout,'(a)') 'Error opening PDB file.'
602 c print *,'Begin reading pdb data'
604 c print *,'Finished reading pdb data'
605 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
607 c itype_pdb(i)=itype(i)
610 c write (iout,'(a,i3)') 'nsup=',nsup
612 c if (nsup.le.(nct-nnt+1)) then
613 c do i=0,nct-nnt+1-nsup
614 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
620 c & 'Error - sequences to be superposed do not match.'
623 c do i=0,nsup-(nct-nnt+1)
624 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
626 c nstart_sup=nstart_sup+i
632 c & 'Error - sequences to be superposed do not match.'
635 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
636 c & ' nstart_seq=',nstart_seq
641 C write (iout,'(/a,i3,a)')
642 C 'The chain contains',ns,' disulfide-bridging cysteines.'
643 write (iout,'(20i4)') (iss(i),i=1,ns)
645 write(iout,*)"Running with dynamic disulfide-bond formation"
647 write (iout,'(/a/)') 'Pre-formed links are:'
653 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
654 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
660 if (ns.gt.0.and.dyn_ss) then
664 forcon(i-nss)=forcon(i)
671 dyn_ss_mask(iss(i))=.true.
674 c Read distance restraints
675 if (constr_dist.gt.0) then
676 call read_dist_constr
681 c-----------------------------------------------------------------------------
682 logical function seq_comp(itypea,itypeb,length)
684 integer length,itypea(length),itypeb(length)
687 if (itypea(i).ne.itypeb(i)) then
695 c-----------------------------------------------------------------------------
696 subroutine read_bridge
697 C Read information about disulfide bridges.
700 include 'COMMON.IOUNITS'
703 include 'COMMON.INTERACT'
704 include 'COMMON.LOCAL'
705 include 'COMMON.NAMES'
706 include 'COMMON.CHAIN'
707 include 'COMMON.FFIELD'
708 include 'COMMON.SBRIDGE'
709 include 'COMMON.HEADER'
710 include 'COMMON.CONTROL'
711 include 'COMMON.TIME1'
713 include 'COMMON.INFO'
716 C Read bridging residues.
717 read (inp,*) ns,(iss(i),i=1,ns)
719 C Check whether the specified bridging residues are cystines.
721 if (itype(iss(i)).ne.1) then
722 write (iout,'(2a,i3,a)')
723 & 'Do you REALLY think that the residue ',
724 & restyp(itype(iss(i))),i,
725 & ' can form a disulfide bridge?!!!'
726 write (*,'(2a,i3,a)')
727 & 'Do you REALLY think that the residue ',
728 & restyp(itype(iss(i))),i,
729 & ' can form a disulfide bridge?!!!'
731 call mp_stopall(error_msg)
737 C Read preformed bridges.
739 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
742 C Check if the residues involved in bridges are in the specified list of
746 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
747 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
748 write (iout,'(a,i3,a)') 'Disulfide pair',i,
749 & ' contains residues present in other pairs.'
750 write (*,'(a,i3,a)') 'Disulfide pair',i,
751 & ' contains residues present in other pairs.'
753 call mp_stopall(error_msg)
760 if (ihpb(i).eq.iss(j)) goto 10
762 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
765 if (jhpb(i).eq.iss(j)) goto 20
767 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
780 c----------------------------------------------------------------------------
781 subroutine read_angles(kanal,*)
786 include 'COMMON.CHAIN'
787 include 'COMMON.IOUNITS'
789 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
790 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
791 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
792 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
794 theta(i)=deg2rad*theta(i)
795 phi(i)=deg2rad*phi(i)
796 alph(i)=deg2rad*alph(i)
797 omeg(i)=deg2rad*omeg(i)
802 c----------------------------------------------------------------------------
803 subroutine reada(rekord,lancuch,wartosc,default)
805 character*(*) rekord,lancuch
806 double precision wartosc,default
809 iread=index(rekord,lancuch)
814 iread=iread+ilen(lancuch)+1
815 read (rekord(iread:),*) wartosc
818 c----------------------------------------------------------------------------
819 subroutine multreada(rekord,lancuch,tablica,dim,default)
822 double precision tablica(dim),default
823 character*(*) rekord,lancuch
829 iread=index(rekord,lancuch)
830 if (iread.eq.0) return
831 iread=iread+ilen(lancuch)+1
832 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
835 c----------------------------------------------------------------------------
836 subroutine readi(rekord,lancuch,wartosc,default)
838 character*(*) rekord,lancuch
839 integer wartosc,default
842 iread=index(rekord,lancuch)
847 iread=iread+ilen(lancuch)+1
848 read (rekord(iread:),*) wartosc
851 C----------------------------------------------------------------------
852 subroutine multreadi(rekord,lancuch,tablica,dim,default)
855 integer tablica(dim),default
856 character*(*) rekord,lancuch
863 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
864 if (iread.eq.0) return
865 iread=iread+ilen(lancuch)+1
866 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
870 c----------------------------------------------------------------------------
871 subroutine card_concat(card)
873 include 'COMMON.IOUNITS'
875 character*80 karta,ucase
877 read (inp,'(a)') karta
880 do while (karta(80:80).eq.'&')
881 card=card(:ilen(card)+1)//karta(:79)
882 read (inp,'(a)') karta
885 card=card(:ilen(card)+1)//karta
888 c----------------------------------------------------------------------------
897 include 'COMMON.IOUNITS'
898 include 'COMMON.CONTROL'
899 integer lenpre,lenpot,ilen
901 character*16 cformat,cprint
903 integer lenint,lenout
904 call getenv('INPUT',prefix)
905 call getenv('OUTPUT',prefout)
906 call getenv('INTIN',prefintin)
907 call getenv('COORD',cformat)
908 call getenv('PRINTCOOR',cprint)
909 call getenv('SCRATCHDIR',scratchdir)
912 if (index(ucase(cformat),'CX').gt.0) then
919 lenint=ilen(prefintin)
920 C Get the names and open the input files
921 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
923 write (liczba,'(bz,i3.3)') me
924 outname=prefout(:lenout)//'_clust.out_'//liczba
926 outname=prefout(:lenout)//'_clust.out'
929 intinname=prefintin(:lenint)//'.bx'
930 else if (from_cx) then
931 intinname=prefintin(:lenint)//'.cx'
933 intinname=prefintin(:lenint)//'.int'
935 rmsname=prefintin(:lenint)//'.rms'
936 open (jplot,file=prefout(:ilen(prefout))//'.tex',
938 open (jrms,file=rmsname,status='unknown')
939 open(iout,file=outname,status='unknown')
940 C Get parameter filenames and open the parameter files.
941 call getenv('BONDPAR',bondname)
942 open (ibond,file=bondname,status='old')
943 call getenv('THETPAR',thetname)
944 open (ithep,file=thetname,status='old')
945 call getenv('ROTPAR',rotname)
946 open (irotam,file=rotname,status='old')
947 call getenv('TORPAR',torname)
948 open (itorp,file=torname,status='old')
950 call getenv('TORDPAR',tordname)
951 open (itordp,file=tordname,status='old')
953 call getenv('FOURIER',fouriername)
954 open (ifourier,file=fouriername,status='old')
955 call getenv('ELEPAR',elename)
956 open (ielep,file=elename,status='old')
957 call getenv('SIDEPAR',sidename)
958 open (isidep,file=sidename,status='old')
959 call getenv('SIDEP',sidepname)
960 open (isidep1,file=sidepname,status="old")
961 call getenv('SCCORPAR',sccorname)
962 open (isccor,file=sccorname,status="old")
963 call getenv('LIPTRANPAR',liptranname)
964 open (iliptranpar,file=liptranname,status='old')
967 C 8/9/01 In the newest version SCp interaction constants are read from a file
968 C Use -DOLDSCP to use hard-coded constants instead.
970 call getenv('SCPPAR',scpname)
971 open (iscpp,file=scpname,status='old')
975 c--------------------------------------------------------------------------
976 subroutine read_dist_constr
977 implicit real*8 (a-h,o-z)
982 include 'COMMON.CONTROL'
983 include 'COMMON.CHAIN'
984 include 'COMMON.IOUNITS'
985 include 'COMMON.SBRIDGE'
986 include 'COMMON.INTERACT'
987 integer ifrag_(2,100),ipair_(2,100)
988 double precision wfrag_(100),wpair_(100)
989 character*500 controlcard
990 logical lprn /.true./
991 logical normalize,next
993 double precision scal_bfac
994 double precision xlink(4,0:4) /
996 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
997 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
998 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
999 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
1000 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
1001 write (iout,*) "Calling read_dist_constr"
1002 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1008 call card_concat(controlcard)
1009 next = index(controlcard,"NEXT").gt.0
1010 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
1011 write (iout,*) "restr_type",restr_type
1012 call readi(controlcard,"NFRAG",nfrag_,0)
1013 call readi(controlcard,"NFRAG",nfrag_,0)
1014 call readi(controlcard,"NPAIR",npair_,0)
1015 call readi(controlcard,"NDIST",ndist_,0)
1016 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1017 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
1018 if (restr_type.eq.10)
1019 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
1020 if (restr_type.eq.12)
1021 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
1022 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1023 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1024 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1025 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1026 normalize = index(controlcard,"NORMALIZE").gt.0
1027 write (iout,*) "WBOLTZD",wboltzd
1028 write (iout,*) "SCAL_PEAK",scal_peak
1029 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1030 write (iout,*) "IFRAG"
1032 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1034 write (iout,*) "IPAIR"
1036 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1038 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
1040 read(inp,'(a)') pdbfile
1042 & "Distance restraints will be constructed from structure ",pdbfile
1043 open(ipdbin,file=pdbfile,status='old',err=11)
1044 call readpdb(.true.)
1049 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1050 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1051 & ifrag_(2,i)=nstart_sup+nsup-1
1052 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1054 if (wfrag_(i).eq.0.0d0) cycle
1055 do j=ifrag_(1,i),ifrag_(2,i)-1
1056 do k=j+1,ifrag_(2,i)
1057 c write (iout,*) "j",j," k",k
1059 if (restr_type.eq.1) then
1065 forcon(nhpb)=wfrag_(i)
1066 else if (constr_dist.eq.2) then
1067 if (ddjk.le.dist_cut) then
1073 forcon(nhpb)=wfrag_(i)
1075 else if (restr_type.eq.3) then
1081 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1083 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
1084 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1089 if (wpair_(i).eq.0.0d0) cycle
1097 do j=ifrag_(1,ii),ifrag_(2,ii)
1098 do k=ifrag_(1,jj),ifrag_(2,jj)
1100 if (restr_type.eq.1) then
1106 forcon(nhpb)=wpair_(i)
1107 else if (constr_dist.eq.2) then
1108 if (ddjk.le.dist_cut) then
1114 forcon(nhpb)=wpair_(i)
1116 else if (restr_type.eq.3) then
1122 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1124 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1125 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1131 write (iout,*) "Distance restraints as read from input"
1133 if (restr_type.eq.12) then
1134 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1135 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1136 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1137 & fordepth_peak(nhpb_peak+1),npeak
1138 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1139 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1140 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1141 c & fordepth_peak(nhpb_peak+1),npeak
1142 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1143 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1144 nhpb_peak=nhpb_peak+1
1145 irestr_type_peak(nhpb_peak)=12
1146 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1148 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1149 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1150 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1151 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1152 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1153 if (ibecarb_peak(nhpb_peak).eq.3) then
1154 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1155 else if (ibecarb_peak(nhpb_peak).eq.2) then
1156 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1157 else if (ibecarb_peak(nhpb_peak).eq.1) then
1158 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1159 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1161 else if (restr_type.eq.11) then
1162 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1163 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1164 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1165 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1167 irestr_type(nhpb)=11
1168 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1169 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1170 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1171 c if (ibecarb(nhpb).gt.0) then
1172 c ihpb(nhpb)=ihpb(nhpb)+nres
1173 c jhpb(nhpb)=jhpb(nhpb)+nres
1175 if (ibecarb(nhpb).eq.3) then
1176 ihpb(nhpb)=ihpb(nhpb)+nres
1177 else if (ibecarb(nhpb).eq.2) then
1178 ihpb(nhpb)=ihpb(nhpb)+nres
1179 else if (ibecarb(nhpb).eq.1) then
1180 ihpb(nhpb)=ihpb(nhpb)+nres
1181 jhpb(nhpb)=jhpb(nhpb)+nres
1183 else if (restr_type.eq.10) then
1184 c Cross-lonk Markov-like potential
1185 call card_concat(controlcard)
1186 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1187 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1189 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1190 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1191 if (index(controlcard,"ZL").gt.0) then
1193 else if (index(controlcard,"ADH").gt.0) then
1195 else if (index(controlcard,"PDH").gt.0) then
1197 else if (index(controlcard,"DSS").gt.0) then
1202 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1203 & xlink(1,link_type))
1204 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1205 & xlink(2,link_type))
1206 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1207 & xlink(3,link_type))
1208 call reada(controlcard,"SIGMA",forcon(nhpb+1),
1209 & xlink(4,link_type))
1210 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1211 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1212 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1213 if (forcon(nhpb+1).le.0.0d0 .or.
1214 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1216 irestr_type(nhpb)=10
1217 c if (ibecarb(nhpb).gt.0) then
1218 c ihpb(nhpb)=ihpb(nhpb)+nres
1219 c jhpb(nhpb)=jhpb(nhpb)+nres
1221 if (ibecarb(nhpb).eq.3) then
1222 jhpb(nhpb)=jhpb(nhpb)+nres
1223 else if (ibecarb(nhpb).eq.2) then
1224 ihpb(nhpb)=ihpb(nhpb)+nres
1225 else if (ibecarb(nhpb).eq.1) then
1226 ihpb(nhpb)=ihpb(nhpb)+nres
1227 jhpb(nhpb)=jhpb(nhpb)+nres
1229 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1230 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1231 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1235 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1236 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1237 if (forcon(nhpb+1).gt.0.0d0) then
1239 if (dhpb1(nhpb).eq.0.0d0) then
1244 c if (ibecarb(nhpb).gt.0) then
1245 c ihpb(nhpb)=ihpb(nhpb)+nres
1246 c jhpb(nhpb)=jhpb(nhpb)+nres
1248 if (ibecarb(nhpb).eq.3) then
1249 jhpb(nhpb)=jhpb(nhpb)+nres
1250 else if (ibecarb(nhpb).eq.2) then
1251 ihpb(nhpb)=ihpb(nhpb)+nres
1252 else if (ibecarb(nhpb).eq.1) then
1253 ihpb(nhpb)=ihpb(nhpb)+nres
1254 jhpb(nhpb)=jhpb(nhpb)+nres
1256 if (dhpb(nhpb).eq.0.0d0)
1257 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1259 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1260 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1262 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1263 C if (forcon(nhpb+1).gt.0.0d0) then
1265 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1268 if (restr_type.eq.4) then
1269 write (iout,*) "The BFAC array"
1271 write (iout,'(i5,f10.5)') i,bfac(i)
1274 if (itype(i).eq.ntyp1) cycle
1276 if (itype(j).eq.ntyp1) cycle
1277 if (itype(i).eq.10) then
1282 if (itype(j).eq.10) then
1292 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1295 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1296 ihpb(nhpb)=i+nres*ii
1297 jhpb(nhpb)=j+nres*jj
1298 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1299 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1300 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1301 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1310 if (restr_type.eq.5) then
1311 restr_on_coord=.true.
1313 if (itype(i).eq.ntyp1) cycle
1314 bfac(i)=(scal_bfac/bfac(i))**2
1323 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1324 & fordepthmax=fordepth(i)
1327 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1330 if (nhpb.gt.nss) then
1331 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1332 & "The following",nhpb-nss,
1333 & " distance restraints have been imposed:",
1334 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1337 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1338 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1345 11 write (iout,*)"read_dist_restr: error reading reference structure"
1348 c-------------------------------------------------------------------------------
1349 subroutine read_saxs_constr
1350 implicit real*8 (a-h,o-z)
1351 include 'DIMENSIONS'
1355 include 'COMMON.CONTROL'
1356 include 'COMMON.CHAIN'
1357 include 'COMMON.IOUNITS'
1358 include 'COMMON.SBRIDGE'
1359 include 'COMMON.SAXS'
1360 double precision cm(3)
1362 write (iout,*) "Calling read_saxs nsaxs",nsaxs
1364 if (saxs_mode.eq.0) then
1365 c SAXS distance distribution
1367 read(inp,*) distsaxs(i),Psaxs(i)
1371 Cnorm = Cnorm + Psaxs(i)
1373 write (iout,*) "Cnorm",Cnorm
1375 Psaxs(i)=Psaxs(i)/Cnorm
1377 write (iout,*) "Normalized distance distribution from SAXS"
1379 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1383 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1385 write (iout,*) "Wsaxs0",Wsaxs0
1389 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1396 cm(j)=cm(j)+Csaxs(j,i)
1404 Csaxs(j,i)=Csaxs(j,i)-cm(j)
1407 write (iout,*) "SAXS sphere coordinates"
1409 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)