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 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
725 C Check whether the specified bridging residues are cystines.
727 if (itype(iss(i)).ne.1) then
728 write (iout,'(2a,i3,a)')
729 & 'Do you REALLY think that the residue ',
730 & restyp(itype(iss(i))),i,
731 & ' can form a disulfide bridge?!!!'
732 write (*,'(2a,i3,a)')
733 & 'Do you REALLY think that the residue ',
734 & restyp(itype(iss(i))),i,
735 & ' can form a disulfide bridge?!!!'
737 call mp_stopall(error_msg)
743 C Read preformed bridges.
745 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
748 C Check if the residues involved in bridges are in the specified list of
752 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
753 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
754 write (iout,'(a,i3,a)') 'Disulfide pair',i,
755 & ' contains residues present in other pairs.'
756 write (*,'(a,i3,a)') 'Disulfide pair',i,
757 & ' contains residues present in other pairs.'
759 call mp_stopall(error_msg)
766 if (ihpb(i).eq.iss(j)) goto 10
768 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
771 if (jhpb(i).eq.iss(j)) goto 20
773 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
786 c----------------------------------------------------------------------------
787 subroutine read_angles(kanal,*)
792 include 'COMMON.CHAIN'
793 include 'COMMON.IOUNITS'
795 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
796 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
797 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
798 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
800 theta(i)=deg2rad*theta(i)
801 phi(i)=deg2rad*phi(i)
802 alph(i)=deg2rad*alph(i)
803 omeg(i)=deg2rad*omeg(i)
808 c----------------------------------------------------------------------------
809 subroutine reada(rekord,lancuch,wartosc,default)
811 character*(*) rekord,lancuch
812 double precision wartosc,default
815 iread=index(rekord,lancuch)
820 iread=iread+ilen(lancuch)+1
821 read (rekord(iread:),*) wartosc
824 c----------------------------------------------------------------------------
825 subroutine multreada(rekord,lancuch,tablica,dim,default)
828 double precision tablica(dim),default
829 character*(*) rekord,lancuch
835 iread=index(rekord,lancuch)
836 if (iread.eq.0) return
837 iread=iread+ilen(lancuch)+1
838 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
841 c----------------------------------------------------------------------------
842 subroutine readi(rekord,lancuch,wartosc,default)
844 character*(*) rekord,lancuch
845 integer wartosc,default
848 iread=index(rekord,lancuch)
853 iread=iread+ilen(lancuch)+1
854 read (rekord(iread:),*) wartosc
857 C----------------------------------------------------------------------
858 subroutine multreadi(rekord,lancuch,tablica,dim,default)
861 integer tablica(dim),default
862 character*(*) rekord,lancuch
869 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
870 if (iread.eq.0) return
871 iread=iread+ilen(lancuch)+1
872 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
876 c----------------------------------------------------------------------------
877 subroutine card_concat(card)
879 include 'COMMON.IOUNITS'
881 character*80 karta,ucase
883 read (inp,'(a)') karta
886 do while (karta(80:80).eq.'&')
887 card=card(:ilen(card)+1)//karta(:79)
888 read (inp,'(a)') karta
891 card=card(:ilen(card)+1)//karta
894 c----------------------------------------------------------------------------
903 include 'COMMON.IOUNITS'
904 include 'COMMON.CONTROL'
905 integer lenpre,lenpot,ilen
907 character*16 cformat,cprint
909 integer lenint,lenout
910 call getenv('INPUT',prefix)
911 call getenv('OUTPUT',prefout)
912 call getenv('INTIN',prefintin)
913 call getenv('COORD',cformat)
914 call getenv('PRINTCOOR',cprint)
915 call getenv('SCRATCHDIR',scratchdir)
918 if (index(ucase(cformat),'CX').gt.0) then
925 lenint=ilen(prefintin)
926 C Get the names and open the input files
927 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
929 write (liczba,'(bz,i3.3)') me
930 outname=prefout(:lenout)//'_clust.out_'//liczba
932 outname=prefout(:lenout)//'_clust.out'
935 intinname=prefintin(:lenint)//'.bx'
936 else if (from_cx) then
937 intinname=prefintin(:lenint)//'.cx'
939 intinname=prefintin(:lenint)//'.int'
941 rmsname=prefintin(:lenint)//'.rms'
942 open (jplot,file=prefout(:ilen(prefout))//'.tex',
944 open (jrms,file=rmsname,status='unknown')
945 open(iout,file=outname,status='unknown')
946 C Get parameter filenames and open the parameter files.
947 call getenv('BONDPAR',bondname)
948 open (ibond,file=bondname,status='old')
949 call getenv('THETPAR',thetname)
950 open (ithep,file=thetname,status='old')
951 call getenv('ROTPAR',rotname)
952 open (irotam,file=rotname,status='old')
953 call getenv('TORPAR',torname)
954 open (itorp,file=torname,status='old')
956 call getenv('TORDPAR',tordname)
957 open (itordp,file=tordname,status='old')
959 call getenv('FOURIER',fouriername)
960 open (ifourier,file=fouriername,status='old')
961 call getenv('ELEPAR',elename)
962 open (ielep,file=elename,status='old')
963 call getenv('SIDEPAR',sidename)
964 open (isidep,file=sidename,status='old')
965 call getenv('SIDEP',sidepname)
966 open (isidep1,file=sidepname,status="old")
967 call getenv('SCCORPAR',sccorname)
968 open (isccor,file=sccorname,status="old")
969 call getenv('LIPTRANPAR',liptranname)
970 open (iliptranpar,file=liptranname,status='old')
973 C 8/9/01 In the newest version SCp interaction constants are read from a file
974 C Use -DOLDSCP to use hard-coded constants instead.
976 call getenv('SCPPAR',scpname)
977 open (iscpp,file=scpname,status='old')
981 c--------------------------------------------------------------------------
982 subroutine read_dist_constr
983 implicit real*8 (a-h,o-z)
988 include 'COMMON.CONTROL'
989 include 'COMMON.CHAIN'
990 include 'COMMON.IOUNITS'
991 include 'COMMON.SBRIDGE'
992 include 'COMMON.INTERACT'
993 integer ifrag_(2,100),ipair_(2,100)
994 double precision wfrag_(100),wpair_(100)
995 character*500 controlcard
996 logical lprn /.true./
997 logical normalize,next
999 double precision scal_bfac
1000 double precision xlink(4,0:4) /
1002 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
1003 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
1004 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
1005 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
1006 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
1007 write (iout,*) "Calling read_dist_constr"
1008 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1014 call card_concat(controlcard)
1015 next = index(controlcard,"NEXT").gt.0
1016 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
1017 write (iout,*) "restr_type",restr_type
1018 call readi(controlcard,"NFRAG",nfrag_,0)
1019 call readi(controlcard,"NFRAG",nfrag_,0)
1020 call readi(controlcard,"NPAIR",npair_,0)
1021 call readi(controlcard,"NDIST",ndist_,0)
1022 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1023 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
1024 if (restr_type.eq.10)
1025 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
1026 if (restr_type.eq.12)
1027 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
1028 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1029 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1030 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1031 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1032 normalize = index(controlcard,"NORMALIZE").gt.0
1033 write (iout,*) "WBOLTZD",wboltzd
1034 write (iout,*) "SCAL_PEAK",scal_peak
1035 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1036 write (iout,*) "IFRAG"
1038 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1040 write (iout,*) "IPAIR"
1042 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1044 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
1046 read(inp,'(a)') pdbfile
1048 & "Distance restraints will be constructed from structure ",pdbfile
1049 open(ipdbin,file=pdbfile,status='old',err=11)
1050 call readpdb(.true.)
1055 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1056 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1057 & ifrag_(2,i)=nstart_sup+nsup-1
1058 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1060 if (wfrag_(i).eq.0.0d0) cycle
1061 do j=ifrag_(1,i),ifrag_(2,i)-1
1062 do k=j+1,ifrag_(2,i)
1063 c write (iout,*) "j",j," k",k
1065 if (restr_type.eq.1) then
1071 forcon(nhpb)=wfrag_(i)
1072 else if (constr_dist.eq.2) then
1073 if (ddjk.le.dist_cut) then
1079 forcon(nhpb)=wfrag_(i)
1081 else if (restr_type.eq.3) then
1087 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1089 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
1090 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1095 if (wpair_(i).eq.0.0d0) cycle
1103 do j=ifrag_(1,ii),ifrag_(2,ii)
1104 do k=ifrag_(1,jj),ifrag_(2,jj)
1106 if (restr_type.eq.1) then
1112 forcon(nhpb)=wpair_(i)
1113 else if (constr_dist.eq.2) then
1114 if (ddjk.le.dist_cut) then
1120 forcon(nhpb)=wpair_(i)
1122 else if (restr_type.eq.3) then
1128 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1130 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1131 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1137 write (iout,*) "Distance restraints as read from input"
1139 if (restr_type.eq.12) then
1140 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1141 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1142 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1143 & fordepth_peak(nhpb_peak+1),npeak
1144 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1145 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1146 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1147 c & fordepth_peak(nhpb_peak+1),npeak
1148 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1149 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1150 nhpb_peak=nhpb_peak+1
1151 irestr_type_peak(nhpb_peak)=12
1152 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1154 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1155 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1156 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1157 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1158 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1159 if (ibecarb_peak(nhpb_peak).eq.3) then
1160 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1161 else if (ibecarb_peak(nhpb_peak).eq.2) then
1162 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1163 else if (ibecarb_peak(nhpb_peak).eq.1) then
1164 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1165 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1167 else if (restr_type.eq.11) then
1168 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1169 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1170 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1171 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1173 irestr_type(nhpb)=11
1174 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1175 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1176 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1177 c if (ibecarb(nhpb).gt.0) then
1178 c ihpb(nhpb)=ihpb(nhpb)+nres
1179 c jhpb(nhpb)=jhpb(nhpb)+nres
1181 if (ibecarb(nhpb).eq.3) then
1182 ihpb(nhpb)=ihpb(nhpb)+nres
1183 else if (ibecarb(nhpb).eq.2) then
1184 ihpb(nhpb)=ihpb(nhpb)+nres
1185 else if (ibecarb(nhpb).eq.1) then
1186 ihpb(nhpb)=ihpb(nhpb)+nres
1187 jhpb(nhpb)=jhpb(nhpb)+nres
1189 else if (restr_type.eq.10) then
1190 c Cross-lonk Markov-like potential
1191 call card_concat(controlcard)
1192 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1193 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1195 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1196 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1197 if (index(controlcard,"ZL").gt.0) then
1199 else if (index(controlcard,"ADH").gt.0) then
1201 else if (index(controlcard,"PDH").gt.0) then
1203 else if (index(controlcard,"DSS").gt.0) then
1208 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1209 & xlink(1,link_type))
1210 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1211 & xlink(2,link_type))
1212 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1213 & xlink(3,link_type))
1214 call reada(controlcard,"SIGMA",forcon(nhpb+1),
1215 & xlink(4,link_type))
1216 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1217 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1218 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1219 if (forcon(nhpb+1).le.0.0d0 .or.
1220 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1222 irestr_type(nhpb)=10
1223 c if (ibecarb(nhpb).gt.0) then
1224 c ihpb(nhpb)=ihpb(nhpb)+nres
1225 c jhpb(nhpb)=jhpb(nhpb)+nres
1227 if (ibecarb(nhpb).eq.3) then
1228 jhpb(nhpb)=jhpb(nhpb)+nres
1229 else if (ibecarb(nhpb).eq.2) then
1230 ihpb(nhpb)=ihpb(nhpb)+nres
1231 else if (ibecarb(nhpb).eq.1) then
1232 ihpb(nhpb)=ihpb(nhpb)+nres
1233 jhpb(nhpb)=jhpb(nhpb)+nres
1235 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1236 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1237 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1241 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1242 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1243 if (forcon(nhpb+1).gt.0.0d0) then
1245 if (dhpb1(nhpb).eq.0.0d0) then
1250 c if (ibecarb(nhpb).gt.0) then
1251 c ihpb(nhpb)=ihpb(nhpb)+nres
1252 c jhpb(nhpb)=jhpb(nhpb)+nres
1254 if (ibecarb(nhpb).eq.3) then
1255 jhpb(nhpb)=jhpb(nhpb)+nres
1256 else if (ibecarb(nhpb).eq.2) then
1257 ihpb(nhpb)=ihpb(nhpb)+nres
1258 else if (ibecarb(nhpb).eq.1) then
1259 ihpb(nhpb)=ihpb(nhpb)+nres
1260 jhpb(nhpb)=jhpb(nhpb)+nres
1262 if (dhpb(nhpb).eq.0.0d0)
1263 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1265 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1266 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1268 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1269 C if (forcon(nhpb+1).gt.0.0d0) then
1271 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1274 if (restr_type.eq.4) then
1275 write (iout,*) "The BFAC array"
1277 write (iout,'(i5,f10.5)') i,bfac(i)
1280 if (itype(i).eq.ntyp1) cycle
1282 if (itype(j).eq.ntyp1) cycle
1283 if (itype(i).eq.10) then
1288 if (itype(j).eq.10) then
1298 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1301 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1302 ihpb(nhpb)=i+nres*ii
1303 jhpb(nhpb)=j+nres*jj
1304 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1305 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1306 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1307 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1316 if (restr_type.eq.5) then
1317 restr_on_coord=.true.
1319 if (itype(i).eq.ntyp1) cycle
1320 bfac(i)=(scal_bfac/bfac(i))**2
1329 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1330 & fordepthmax=fordepth(i)
1333 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1336 if (nhpb.gt.nss) then
1337 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1338 & "The following",nhpb-nss,
1339 & " distance restraints have been imposed:",
1340 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1343 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1344 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1351 11 write (iout,*)"read_dist_restr: error reading reference structure"
1354 c-------------------------------------------------------------------------------
1355 subroutine read_saxs_constr
1356 implicit real*8 (a-h,o-z)
1357 include 'DIMENSIONS'
1361 include 'COMMON.CONTROL'
1362 include 'COMMON.CHAIN'
1363 include 'COMMON.IOUNITS'
1364 include 'COMMON.SBRIDGE'
1365 include 'COMMON.SAXS'
1366 double precision cm(3)
1368 write (iout,*) "Calling read_saxs nsaxs",nsaxs
1370 if (saxs_mode.eq.0) then
1371 c SAXS distance distribution
1373 read(inp,*) distsaxs(i),Psaxs(i)
1377 Cnorm = Cnorm + Psaxs(i)
1379 write (iout,*) "Cnorm",Cnorm
1381 Psaxs(i)=Psaxs(i)/Cnorm
1383 write (iout,*) "Normalized distance distribution from SAXS"
1385 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
1389 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
1391 write (iout,*) "Wsaxs0",Wsaxs0
1395 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
1402 cm(j)=cm(j)+Csaxs(j,i)
1410 Csaxs(j,i)=Csaxs(j,i)-cm(j)
1413 write (iout,*) "SAXS sphere coordinates"
1415 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)