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'
20 character*320 controlcard,ucase
24 integer i,i1,i2,it1,it2
26 read (INP,'(a80)') titel
27 call card_concat(controlcard)
29 call readi(controlcard,'NRES',nres,0)
30 call readi(controlcard,'RESCALE',rescale_mode,2)
31 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
32 write (iout,*) "DISTCHAINMAX",distchainmax
33 C Reading the dimensions of box in x,y,z coordinates
34 call reada(controlcard,'BOXX',boxxsize,100.0d0)
35 call reada(controlcard,'BOXY',boxysize,100.0d0)
36 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
37 c Cutoff range for interactions
38 call reada(controlcard,"R_CUT",r_cut,15.0d0)
39 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
40 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
41 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
42 if (lipthick.gt.0.0d0) then
43 bordliptop=(boxzsize+lipthick)/2.0
44 bordlipbot=bordliptop-lipthick
46 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
47 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
48 buflipbot=bordlipbot+lipbufthick
49 bufliptop=bordliptop-lipbufthick
50 if ((lipbufthick*2.0d0).gt.lipthick)
51 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
53 write(iout,*) "bordliptop=",bordliptop
54 write(iout,*) "bordlipbot=",bordlipbot
55 write(iout,*) "bufliptop=",bufliptop
56 write(iout,*) "buflipbot=",buflipbot
58 call readi(controlcard,'SHIELD',shield_mode,0)
59 write (iout,*) "SHIELD MODE",shield_mode
60 call readi(controlcard,'TUBEMOD',tubelog,0)
61 write (iout,*) TUBElog,"TUBEMODE"
62 if (shield_mode.gt.0) then
64 C VSolvSphere the volume of solving sphere
66 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
67 C there will be no distinction between proline peptide group and normal peptide
68 C group in case of shielding parameters
69 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
70 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
71 write (iout,*) VSolvSphere,VSolvSphere_div
72 C long axis of side chain
74 long_r_sidechain(i)=vbldsc0(1,i)
75 short_r_sidechain(i)=sigma0(i)
79 if (TUBElog.gt.0) then
80 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
81 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
82 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
83 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
84 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
85 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
86 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
87 buftubebot=bordtubebot+tubebufthick
88 buftubetop=bordtubetop-tubebufthick
90 call readi(controlcard,'PDBOUT',outpdb,0)
91 call readi(controlcard,'MOL2OUT',outmol2,0)
92 refstr=(index(controlcard,'REFSTR').gt.0)
93 write (iout,*) "REFSTR",refstr
94 pdbref=(index(controlcard,'PDBREF').gt.0)
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 call readi(controlcard,'NCUT',ncut,1)
108 call readi(controlcard,'SYM',symetr,1)
109 write (iout,*) 'sym', symetr
110 call readi(controlcard,'NSTART',nstart,0)
111 call readi(controlcard,'NEND',nend,0)
112 call reada(controlcard,'ECUT',ecut,10.0d0)
113 call reada(controlcard,'PROB',prob_limit,0.99d0)
114 write (iout,*) "Probability limit",prob_limit
115 lgrp=(index(controlcard,'LGRP').gt.0)
116 caonly=(index(controlcard,'CA_ONLY').gt.0)
117 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
118 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
119 call readi(controlcard,'IOPT',iopt,2)
120 lside = index(controlcard,"SIDE").gt.0
121 efree = index(controlcard,"EFREE").gt.0
122 call readi(controlcard,'NTEMP',nT,1)
123 write (iout,*) "nT",nT
124 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
125 write (iout,*) "nT",nT
126 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
128 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
130 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
131 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
132 lprint_int=index(controlcard,"PRINT_INT") .gt.0
136 c--------------------------------------------------------------------------
139 C Read molecular data.
143 include 'COMMON.IOUNITS'
146 include 'COMMON.INTERACT'
147 include 'COMMON.LOCAL'
148 include 'COMMON.NAMES'
149 include 'COMMON.CHAIN'
150 include 'COMMON.FFIELD'
151 include 'COMMON.SBRIDGE'
152 include 'COMMON.HEADER'
153 include 'COMMON.CONTROL'
154 include 'COMMON.CONTACTS'
155 include 'COMMON.TIME1'
156 include 'COMMON.TORCNSTR'
157 include 'COMMON.SHIELD'
159 include 'COMMON.INFO'
161 character*4 sequence(maxres)
162 character*800 weightcard
164 double precision x(maxvar)
165 integer itype_pdb(maxres)
167 integer i,j,kkk,i1,i2,it1,it2
171 C Read weights of the subsequent energy terms.
172 call card_concat(weightcard)
173 write(iout,*) weightcard
174 C call reada(weightcard,'WSC',wsc,1.0d0)
176 call reada(weightcard,'WLONG',wsc,wsc)
177 call reada(weightcard,'WSCP',wscp,1.0d0)
178 call reada(weightcard,'WELEC',welec,1.0D0)
179 call reada(weightcard,'WVDWPP',wvdwpp,welec)
180 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
181 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
182 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
183 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
184 call reada(weightcard,'WTURN3',wturn3,1.0D0)
185 call reada(weightcard,'WTURN4',wturn4,1.0D0)
186 call reada(weightcard,'WTURN6',wturn6,1.0D0)
187 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
188 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
189 call reada(weightcard,'WBOND',wbond,1.0D0)
190 call reada(weightcard,'WTOR',wtor,1.0D0)
191 call reada(weightcard,'WTORD',wtor_d,1.0D0)
192 call reada(weightcard,'WANG',wang,1.0D0)
193 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
194 call reada(weightcard,'SCAL14',scal14,0.4D0)
195 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
196 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
197 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
198 if (index(weightcard,'SOFT').gt.0) ipot=6
199 call reada(weightcard,"D0CM",d0cm,3.78d0)
200 call reada(weightcard,"AKCM",akcm,15.1d0)
201 call reada(weightcard,"AKTH",akth,11.0d0)
202 call reada(weightcard,"AKCT",akct,12.0d0)
203 call reada(weightcard,"V1SS",v1ss,-1.08d0)
204 call reada(weightcard,"V2SS",v2ss,7.61d0)
205 call reada(weightcard,"V3SS",v3ss,13.7d0)
206 call reada(weightcard,"EBR",ebr,-5.50D0)
207 call reada(weightcard,'WSHIELD',wshield,1.0d0)
208 write(iout,*) 'WSHIELD',wshield
209 call reada(weightcard,'WTUBE',wtube,0.0d0)
210 write(iout,*) 'WTUBE',wtube
211 call reada(weightcard,'WLT',wliptran,0.0D0)
212 call reada(weightcard,"ATRISS",atriss,0.301D0)
213 call reada(weightcard,"BTRISS",btriss,0.021D0)
214 call reada(weightcard,"CTRISS",ctriss,1.001D0)
215 call reada(weightcard,"DTRISS",dtriss,1.001D0)
216 write (iout,*) "ATRISS=", atriss
217 write (iout,*) "BTRISS=", btriss
218 write (iout,*) "CTRISS=", ctriss
219 write (iout,*) "DTRISS=", dtriss
220 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
222 dyn_ss_mask(i)=.false.
226 dyn_ssbond_ij(i,j)=1.0d300
229 call reada(weightcard,"HT",Ht,0.0D0)
231 ss_depth=ebr/wsc-0.25*eps(1,1)
232 Ht=Ht/wsc-0.25*eps(1,1)
233 akcm=akcm*wstrain/wsc
234 akth=akth*wstrain/wsc
235 akct=akct*wstrain/wsc
236 v1ss=v1ss*wstrain/wsc
237 v2ss=v2ss*wstrain/wsc
238 v3ss=v3ss*wstrain/wsc
240 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
242 write (iout,'(/a)') "Disulfide bridge parameters:"
243 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
244 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
245 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
246 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
249 C 12/1/95 Added weight for the multi-body term WCORR
250 call reada(weightcard,'WCORRH',wcorr,1.0D0)
251 if (wcorr4.gt.0.0d0) wcorr=wcorr4
271 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
272 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
273 & wturn4,wturn6,wsccor
274 10 format (/'Energy-term weights (unscaled):'//
275 & 'WSCC= ',f10.6,' (SC-SC)'/
276 & 'WSCP= ',f10.6,' (SC-p)'/
277 & 'WELEC= ',f10.6,' (p-p electr)'/
278 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
279 & 'WBOND= ',f10.6,' (stretching)'/
280 & 'WANG= ',f10.6,' (bending)'/
281 & 'WSCLOC= ',f10.6,' (SC local)'/
282 & 'WTOR= ',f10.6,' (torsional)'/
283 & 'WTORD= ',f10.6,' (double torsional)'/
284 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
285 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
286 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
287 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
288 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
289 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
290 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
291 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
292 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
294 if (wcorr4.gt.0.0d0) then
295 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
296 & 'between contact pairs of peptide groups'
297 write (iout,'(2(a,f5.3/))')
298 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
299 & 'Range of quenching the correlation terms:',2*delt_corr
300 else if (wcorr.gt.0.0d0) then
301 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
302 & 'between contact pairs of peptide groups'
304 write (iout,'(a,f8.3)')
305 & 'Scaling factor of 1,4 SC-p interactions:',scal14
306 write (iout,'(a,f8.3)')
307 & 'General scaling factor of SC-p interactions:',scalscp
308 r0_corr=cutoff_corr-delt_corr
310 aad(i,1)=scalscp*aad(i,1)
311 aad(i,2)=scalscp*aad(i,2)
312 bad(i,1)=scalscp*bad(i,1)
313 bad(i,2)=scalscp*bad(i,2)
317 print *,'indpdb=',indpdb,' pdbref=',pdbref
319 C Read sequence if not taken from the pdb file.
320 if (iscode.gt.0) then
321 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
323 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
325 C Convert sequence to numeric code
327 itype(i)=rescode(i,sequence(i),iscode)
330 print '(20i4)',(itype(i),i=1,nres)
334 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
336 if (itype(i).eq.ntyp1) then
340 else if (iabs(itype(i+1)).ne.20) then
342 else if (iabs(itype(i)).ne.20) then
349 write (iout,*) "ITEL"
351 write (iout,*) i,itype(i),itel(i)
354 print *,'Call Read_Bridge.'
356 C this fragment reads diheadral constrains
357 if (with_dihed_constr) then
359 read (inp,*) ndih_constr
360 if (ndih_constr.gt.0) then
362 C write (iout,*) 'FTORS',ftors
363 C ftors is the force constant for torsional quartic constrains
364 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
367 & 'There are',ndih_constr,' constraints on phi angles.'
369 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
373 phi0(i)=deg2rad*phi0(i)
374 drange(i)=deg2rad*drange(i)
376 endif ! endif ndif_constr.gt.0
377 endif ! with_dihed_constr
378 if (with_theta_constr) then
379 C with_theta_constr is keyword allowing for occurance of theta constrains
380 read (inp,*) ntheta_constr
381 C ntheta_constr is the number of theta constrains
382 if (ntheta_constr.gt.0) then
384 read (inp,*) (itheta_constr(i),theta_constr0(i),
385 & theta_drange(i),for_thet_constr(i),
387 C the above code reads from 1 to ntheta_constr
388 C itheta_constr(i) residue i for which is theta_constr
389 C theta_constr0 the global minimum value
390 C theta_drange is range for which there is no energy penalty
391 C for_thet_constr is the force constant for quartic energy penalty
393 C if(me.eq.king.or..not.out1file)then
395 & 'There are',ntheta_constr,' constraints on phi angles.'
397 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
403 theta_constr0(i)=deg2rad*theta_constr0(i)
404 theta_drange(i)=deg2rad*theta_drange(i)
406 C if(me.eq.king.or..not.out1file)
407 C & write (iout,*) 'FTORS',ftors
408 C do i=1,ntheta_constr
409 C ii = itheta_constr(i)
410 C thetabound(1,ii) = phi0(i)-drange(i)
411 C thetabound(2,ii) = phi0(i)+drange(i)
413 endif ! ntheta_constr.gt.0
414 endif! with_theta_constr
418 print *,'NNT=',NNT,' NCT=',NCT
419 if (itype(1).eq.ntyp1) nnt=2
420 if (itype(nres).eq.ntyp1) nct=nct-1
421 if (nstart.lt.nnt) nstart=nnt
422 if (nend.gt.nct .or. nend.eq.0) nend=nct
423 write (iout,*) "nstart",nstart," nend",nend
426 c read(inp,'(a)') pdbfile
427 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
428 c open(ipdbin,file=pdbfile,status='old',err=33)
430 c 33 write (iout,'(a)') 'Error opening PDB file.'
433 c print *,'Begin reading pdb data'
435 c print *,'Finished reading pdb data'
436 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
438 c itype_pdb(i)=itype(i)
441 c write (iout,'(a,i3)') 'nsup=',nsup
443 c if (nsup.le.(nct-nnt+1)) then
444 c do i=0,nct-nnt+1-nsup
445 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
451 c & 'Error - sequences to be superposed do not match.'
454 c do i=0,nsup-(nct-nnt+1)
455 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
457 c nstart_sup=nstart_sup+i
463 c & 'Error - sequences to be superposed do not match.'
466 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
467 c & ' nstart_seq=',nstart_seq
471 write (iout,*) "molread: REFSTR",refstr
473 if (.not.pdbref) then
474 call read_angles(inp,*38)
476 38 write (iout,'(a)') 'Error reading reference structure.'
478 call mp_stopall(Error_Msg)
480 stop 'Error reading reference structure'
493 call contact(.true.,ncont_ref,icont_ref)
496 C write (iout,'(/a,i3,a)')
497 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
498 write (iout,'(20i4)') (iss(i),i=1,ns)
500 write(iout,*)"Running with dynamic disulfide-bond formation"
502 write (iout,'(/a/)') 'Pre-formed links are:'
508 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
509 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
515 if (ns.gt.0.and.dyn_ss) then
519 forcon(i-nss)=forcon(i)
526 dyn_ss_mask(iss(i))=.true.
529 c Read distance restraints
530 if (constr_dist.gt.0) then
531 call read_dist_constr
536 c-----------------------------------------------------------------------------
537 logical function seq_comp(itypea,itypeb,length)
539 integer length,itypea(length),itypeb(length)
542 if (itypea(i).ne.itypeb(i)) then
550 c-----------------------------------------------------------------------------
551 subroutine read_bridge
552 C Read information about disulfide bridges.
555 include 'COMMON.IOUNITS'
558 include 'COMMON.INTERACT'
559 include 'COMMON.LOCAL'
560 include 'COMMON.NAMES'
561 include 'COMMON.CHAIN'
562 include 'COMMON.FFIELD'
563 include 'COMMON.SBRIDGE'
564 include 'COMMON.HEADER'
565 include 'COMMON.CONTROL'
566 include 'COMMON.TIME1'
568 include 'COMMON.INFO'
571 C Read bridging residues.
572 read (inp,*) ns,(iss(i),i=1,ns)
574 C Check whether the specified bridging residues are cystines.
576 if (itype(iss(i)).ne.1) then
577 write (iout,'(2a,i3,a)')
578 & 'Do you REALLY think that the residue ',
579 & restyp(itype(iss(i))),i,
580 & ' can form a disulfide bridge?!!!'
581 write (*,'(2a,i3,a)')
582 & 'Do you REALLY think that the residue ',
583 & restyp(itype(iss(i))),i,
584 & ' can form a disulfide bridge?!!!'
586 call mp_stopall(error_msg)
592 C Read preformed bridges.
594 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
597 C Check if the residues involved in bridges are in the specified list of
601 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
602 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
603 write (iout,'(a,i3,a)') 'Disulfide pair',i,
604 & ' contains residues present in other pairs.'
605 write (*,'(a,i3,a)') 'Disulfide pair',i,
606 & ' contains residues present in other pairs.'
608 call mp_stopall(error_msg)
615 if (ihpb(i).eq.iss(j)) goto 10
617 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
620 if (jhpb(i).eq.iss(j)) goto 20
622 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
635 c----------------------------------------------------------------------------
636 subroutine read_angles(kanal,*)
641 include 'COMMON.CHAIN'
642 include 'COMMON.IOUNITS'
644 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
645 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
646 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
647 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
649 theta(i)=deg2rad*theta(i)
650 phi(i)=deg2rad*phi(i)
651 alph(i)=deg2rad*alph(i)
652 omeg(i)=deg2rad*omeg(i)
657 c----------------------------------------------------------------------------
658 subroutine reada(rekord,lancuch,wartosc,default)
660 character*(*) rekord,lancuch
661 double precision wartosc,default
664 iread=index(rekord,lancuch)
669 iread=iread+ilen(lancuch)+1
670 read (rekord(iread:),*) wartosc
673 c----------------------------------------------------------------------------
674 subroutine multreada(rekord,lancuch,tablica,dim,default)
677 double precision tablica(dim),default
678 character*(*) rekord,lancuch
684 iread=index(rekord,lancuch)
685 if (iread.eq.0) return
686 iread=iread+ilen(lancuch)+1
687 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
690 c----------------------------------------------------------------------------
691 subroutine readi(rekord,lancuch,wartosc,default)
693 character*(*) rekord,lancuch
694 integer wartosc,default
697 iread=index(rekord,lancuch)
702 iread=iread+ilen(lancuch)+1
703 read (rekord(iread:),*) wartosc
706 C----------------------------------------------------------------------
707 subroutine multreadi(rekord,lancuch,tablica,dim,default)
710 integer tablica(dim),default
711 character*(*) rekord,lancuch
718 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
719 if (iread.eq.0) return
720 iread=iread+ilen(lancuch)+1
721 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
725 c----------------------------------------------------------------------------
726 subroutine card_concat(card)
728 include 'COMMON.IOUNITS'
730 character*80 karta,ucase
732 read (inp,'(a)') karta
735 do while (karta(80:80).eq.'&')
736 card=card(:ilen(card)+1)//karta(:79)
737 read (inp,'(a)') karta
740 card=card(:ilen(card)+1)//karta
743 c----------------------------------------------------------------------------
752 include 'COMMON.IOUNITS'
753 include 'COMMON.CONTROL'
754 integer lenpre,lenpot,ilen
756 character*16 cformat,cprint
758 integer lenint,lenout
759 call getenv('INPUT',prefix)
760 call getenv('OUTPUT',prefout)
761 call getenv('INTIN',prefintin)
762 call getenv('COORD',cformat)
763 call getenv('PRINTCOOR',cprint)
764 call getenv('SCRATCHDIR',scratchdir)
767 if (index(ucase(cformat),'CX').gt.0) then
774 lenint=ilen(prefintin)
775 C Get the names and open the input files
776 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
778 write (liczba,'(bz,i3.3)') me
779 outname=prefout(:lenout)//'_clust.out_'//liczba
781 outname=prefout(:lenout)//'_clust.out'
784 intinname=prefintin(:lenint)//'.bx'
785 else if (from_cx) then
786 intinname=prefintin(:lenint)//'.cx'
788 intinname=prefintin(:lenint)//'.int'
790 rmsname=prefintin(:lenint)//'.rms'
791 open (jplot,file=prefout(:ilen(prefout))//'.tex',
793 open (jrms,file=rmsname,status='unknown')
794 open(iout,file=outname,status='unknown')
795 C Get parameter filenames and open the parameter files.
796 call getenv('BONDPAR',bondname)
797 open (ibond,file=bondname,status='old')
798 call getenv('THETPAR',thetname)
799 open (ithep,file=thetname,status='old')
800 call getenv('ROTPAR',rotname)
801 open (irotam,file=rotname,status='old')
802 call getenv('TORPAR',torname)
803 open (itorp,file=torname,status='old')
804 call getenv('TORDPAR',tordname)
805 open (itordp,file=tordname,status='old')
806 call getenv('FOURIER',fouriername)
807 open (ifourier,file=fouriername,status='old')
808 call getenv('ELEPAR',elename)
809 open (ielep,file=elename,status='old')
810 call getenv('SIDEPAR',sidename)
811 open (isidep,file=sidename,status='old')
812 call getenv('SIDEP',sidepname)
813 open (isidep1,file=sidepname,status="old")
814 call getenv('SCCORPAR',sccorname)
815 open (isccor,file=sccorname,status="old")
816 call getenv('LIPTRANPAR',liptranname)
817 open (iliptranpar,file=liptranname,status='old')
818 call getenv('TUBEPAR',tubename)
819 open (itube,file=tubename,status='old')
823 C 8/9/01 In the newest version SCp interaction constants are read from a file
824 C Use -DOLDSCP to use hard-coded constants instead.
826 call getenv('SCPPAR',scpname)
827 open (iscpp,file=scpname,status='old')
831 subroutine read_dist_constr
832 implicit real*8 (a-h,o-z)
837 include 'COMMON.CONTROL'
838 include 'COMMON.CHAIN'
839 include 'COMMON.IOUNITS'
840 include 'COMMON.SBRIDGE'
841 integer ifrag_(2,100),ipair_(2,100)
842 double precision wfrag_(100),wpair_(100)
843 character*500 controlcard
844 logical lprn /.true./
845 write (iout,*) "Calling read_dist_constr"
846 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
848 write(iout,*) "TU sie wywalam?"
849 call card_concat(controlcard)
850 write (iout,*) controlcard
852 call readi(controlcard,"NFRAG",nfrag_,0)
853 call readi(controlcard,"NPAIR",npair_,0)
854 call readi(controlcard,"NDIST",ndist_,0)
855 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
856 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
857 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
858 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
859 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
860 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
861 write (iout,*) "IFRAG"
863 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
865 write (iout,*) "IPAIR"
867 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
871 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
872 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
873 & ifrag_(2,i)=nstart_sup+nsup-1
874 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
876 if (wfrag_(i).gt.0.0d0) then
877 do j=ifrag_(1,i),ifrag_(2,i)-1
879 write (iout,*) "j",j," k",k
881 if (constr_dist.eq.1) then
886 forcon(nhpb)=wfrag_(i)
887 else if (constr_dist.eq.2) then
888 if (ddjk.le.dist_cut) then
893 forcon(nhpb)=wfrag_(i)
900 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
903 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
904 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
910 if (wpair_(i).gt.0.0d0) then
918 do j=ifrag_(1,ii),ifrag_(2,ii)
919 do k=ifrag_(1,jj),ifrag_(2,jj)
923 forcon(nhpb)=wpair_(i)
925 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
926 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
932 if (constr_dist.eq.11) then
933 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
934 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
935 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
936 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
937 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
939 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
941 if (forcon(nhpb+1).gt.0.0d0) then
943 if (ibecarb(i).gt.0) then
947 if (dhpb(nhpb).eq.0.0d0)
948 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
949 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
950 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
951 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)