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 if (shield_mode.gt.0) then
62 C VSolvSphere the volume of solving sphere
64 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
65 C there will be no distinction between proline peptide group and normal peptide
66 C group in case of shielding parameters
67 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
68 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
69 write (iout,*) VSolvSphere,VSolvSphere_div
70 C long axis of side chain
72 long_r_sidechain(i)=vbldsc0(1,i)
73 short_r_sidechain(i)=sigma0(i)
77 call readi(controlcard,'PDBOUT',outpdb,0)
78 call readi(controlcard,'MOL2OUT',outmol2,0)
79 refstr=(index(controlcard,'REFSTR').gt.0)
80 write (iout,*) "REFSTR",refstr
81 pdbref=(index(controlcard,'PDBREF').gt.0)
82 iscode=index(controlcard,'ONE_LETTER')
83 tree=(index(controlcard,'MAKE_TREE').gt.0)
84 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
85 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
86 write (iout,*) "with_dihed_constr ",with_dihed_constr,
87 & " CONSTR_DIST",constr_dist
88 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
89 write (iout,*) "with_theta_constr ",with_theta_constr
91 min_var=(index(controlcard,'MINVAR').gt.0)
92 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
93 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
94 call readi(controlcard,'NCUT',ncut,1)
95 call readi(controlcard,'SYM',symetr,1)
96 write (iout,*) 'sym', symetr
97 call readi(controlcard,'NSTART',nstart,0)
98 call readi(controlcard,'NEND',nend,0)
99 call reada(controlcard,'ECUT',ecut,10.0d0)
100 call reada(controlcard,'PROB',prob_limit,0.99d0)
101 write (iout,*) "Probability limit",prob_limit
102 lgrp=(index(controlcard,'LGRP').gt.0)
103 caonly=(index(controlcard,'CA_ONLY').gt.0)
104 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
105 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
106 call readi(controlcard,'IOPT',iopt,2)
107 lside = index(controlcard,"SIDE").gt.0
108 efree = index(controlcard,"EFREE").gt.0
109 call readi(controlcard,'NTEMP',nT,1)
110 write (iout,*) "nT",nT
111 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
112 write (iout,*) "nT",nT
113 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
115 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
117 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
118 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
119 lprint_int=index(controlcard,"PRINT_INT") .gt.0
123 c--------------------------------------------------------------------------
126 C Read molecular data.
130 include 'COMMON.IOUNITS'
133 include 'COMMON.INTERACT'
134 include 'COMMON.LOCAL'
135 include 'COMMON.NAMES'
136 include 'COMMON.CHAIN'
137 include 'COMMON.FFIELD'
138 include 'COMMON.SBRIDGE'
139 include 'COMMON.HEADER'
140 include 'COMMON.CONTROL'
141 include 'COMMON.CONTACTS'
142 include 'COMMON.TIME1'
143 include 'COMMON.TORCNSTR'
144 include 'COMMON.SHIELD'
146 include 'COMMON.INFO'
148 character*4 sequence(maxres)
149 character*800 weightcard
151 double precision x(maxvar)
152 integer itype_pdb(maxres)
154 integer i,j,kkk,i1,i2,it1,it2
158 C Read weights of the subsequent energy terms.
159 call card_concat(weightcard)
160 write(iout,*) weightcard
161 C call reada(weightcard,'WSC',wsc,1.0d0)
163 call reada(weightcard,'WLONG',wsc,wsc)
164 call reada(weightcard,'WSCP',wscp,1.0d0)
165 call reada(weightcard,'WELEC',welec,1.0D0)
166 call reada(weightcard,'WVDWPP',wvdwpp,welec)
167 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
168 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
169 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
170 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
171 call reada(weightcard,'WTURN3',wturn3,1.0D0)
172 call reada(weightcard,'WTURN4',wturn4,1.0D0)
173 call reada(weightcard,'WTURN6',wturn6,1.0D0)
174 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
175 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
176 call reada(weightcard,'WBOND',wbond,1.0D0)
177 call reada(weightcard,'WTOR',wtor,1.0D0)
178 call reada(weightcard,'WTORD',wtor_d,1.0D0)
179 call reada(weightcard,'WANG',wang,1.0D0)
180 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
181 call reada(weightcard,'SCAL14',scal14,0.4D0)
182 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
183 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
184 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
185 if (index(weightcard,'SOFT').gt.0) ipot=6
186 call reada(weightcard,"D0CM",d0cm,3.78d0)
187 call reada(weightcard,"AKCM",akcm,15.1d0)
188 call reada(weightcard,"AKTH",akth,11.0d0)
189 call reada(weightcard,"AKCT",akct,12.0d0)
190 call reada(weightcard,"V1SS",v1ss,-1.08d0)
191 call reada(weightcard,"V2SS",v2ss,7.61d0)
192 call reada(weightcard,"V3SS",v3ss,13.7d0)
193 call reada(weightcard,"EBR",ebr,-5.50D0)
194 call reada(weightcard,'WSHIELD',wshield,1.0d0)
195 write(iout,*) 'WSHIELD',wshield
196 call reada(weightcard,'WLT',wliptran,0.0D0)
197 call reada(weightcard,"ATRISS",atriss,0.301D0)
198 call reada(weightcard,"BTRISS",btriss,0.021D0)
199 call reada(weightcard,"CTRISS",ctriss,1.001D0)
200 call reada(weightcard,"DTRISS",dtriss,1.001D0)
201 write (iout,*) "ATRISS=", atriss
202 write (iout,*) "BTRISS=", btriss
203 write (iout,*) "CTRISS=", ctriss
204 write (iout,*) "DTRISS=", dtriss
205 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
207 dyn_ss_mask(i)=.false.
211 dyn_ssbond_ij(i,j)=1.0d300
214 call reada(weightcard,"HT",Ht,0.0D0)
216 ss_depth=ebr/wsc-0.25*eps(1,1)
217 Ht=Ht/wsc-0.25*eps(1,1)
218 akcm=akcm*wstrain/wsc
219 akth=akth*wstrain/wsc
220 akct=akct*wstrain/wsc
221 v1ss=v1ss*wstrain/wsc
222 v2ss=v2ss*wstrain/wsc
223 v3ss=v3ss*wstrain/wsc
225 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
227 write (iout,'(/a)') "Disulfide bridge parameters:"
228 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
229 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
230 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
231 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
234 C 12/1/95 Added weight for the multi-body term WCORR
235 call reada(weightcard,'WCORRH',wcorr,1.0D0)
236 if (wcorr4.gt.0.0d0) wcorr=wcorr4
256 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
257 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
258 & wturn4,wturn6,wsccor
259 10 format (/'Energy-term weights (unscaled):'//
260 & 'WSCC= ',f10.6,' (SC-SC)'/
261 & 'WSCP= ',f10.6,' (SC-p)'/
262 & 'WELEC= ',f10.6,' (p-p electr)'/
263 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
264 & 'WBOND= ',f10.6,' (stretching)'/
265 & 'WANG= ',f10.6,' (bending)'/
266 & 'WSCLOC= ',f10.6,' (SC local)'/
267 & 'WTOR= ',f10.6,' (torsional)'/
268 & 'WTORD= ',f10.6,' (double torsional)'/
269 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
270 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
271 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
272 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
273 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
274 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
275 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
276 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
277 & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
279 if (wcorr4.gt.0.0d0) then
280 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
281 & 'between contact pairs of peptide groups'
282 write (iout,'(2(a,f5.3/))')
283 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
284 & 'Range of quenching the correlation terms:',2*delt_corr
285 else if (wcorr.gt.0.0d0) then
286 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
287 & 'between contact pairs of peptide groups'
289 write (iout,'(a,f8.3)')
290 & 'Scaling factor of 1,4 SC-p interactions:',scal14
291 write (iout,'(a,f8.3)')
292 & 'General scaling factor of SC-p interactions:',scalscp
293 r0_corr=cutoff_corr-delt_corr
295 aad(i,1)=scalscp*aad(i,1)
296 aad(i,2)=scalscp*aad(i,2)
297 bad(i,1)=scalscp*bad(i,1)
298 bad(i,2)=scalscp*bad(i,2)
302 print *,'indpdb=',indpdb,' pdbref=',pdbref
304 C Read sequence if not taken from the pdb file.
305 if (iscode.gt.0) then
306 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
308 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
310 C Convert sequence to numeric code
312 itype(i)=rescode(i,sequence(i),iscode)
315 print '(20i4)',(itype(i),i=1,nres)
319 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
321 if (itype(i).eq.ntyp1) then
325 else if (iabs(itype(i+1)).ne.20) then
327 else if (iabs(itype(i)).ne.20) then
334 write (iout,*) "ITEL"
336 write (iout,*) i,itype(i),itel(i)
339 print *,'Call Read_Bridge.'
341 C this fragment reads diheadral constrains
342 if (with_dihed_constr) then
344 read (inp,*) ndih_constr
345 if (ndih_constr.gt.0) then
347 C write (iout,*) 'FTORS',ftors
348 C ftors is the force constant for torsional quartic constrains
349 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
352 & 'There are',ndih_constr,' constraints on phi angles.'
354 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
358 phi0(i)=deg2rad*phi0(i)
359 drange(i)=deg2rad*drange(i)
361 endif ! endif ndif_constr.gt.0
362 endif ! with_dihed_constr
363 if (with_theta_constr) then
364 C with_theta_constr is keyword allowing for occurance of theta constrains
365 read (inp,*) ntheta_constr
366 C ntheta_constr is the number of theta constrains
367 if (ntheta_constr.gt.0) then
369 read (inp,*) (itheta_constr(i),theta_constr0(i),
370 & theta_drange(i),for_thet_constr(i),
372 C the above code reads from 1 to ntheta_constr
373 C itheta_constr(i) residue i for which is theta_constr
374 C theta_constr0 the global minimum value
375 C theta_drange is range for which there is no energy penalty
376 C for_thet_constr is the force constant for quartic energy penalty
378 C if(me.eq.king.or..not.out1file)then
380 & 'There are',ntheta_constr,' constraints on phi angles.'
382 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
388 theta_constr0(i)=deg2rad*theta_constr0(i)
389 theta_drange(i)=deg2rad*theta_drange(i)
391 C if(me.eq.king.or..not.out1file)
392 C & write (iout,*) 'FTORS',ftors
393 C do i=1,ntheta_constr
394 C ii = itheta_constr(i)
395 C thetabound(1,ii) = phi0(i)-drange(i)
396 C thetabound(2,ii) = phi0(i)+drange(i)
398 endif ! ntheta_constr.gt.0
399 endif! with_theta_constr
403 print *,'NNT=',NNT,' NCT=',NCT
404 if (itype(1).eq.ntyp1) nnt=2
405 if (itype(nres).eq.ntyp1) nct=nct-1
406 if (nstart.lt.nnt) nstart=nnt
407 if (nend.gt.nct .or. nend.eq.0) nend=nct
408 write (iout,*) "nstart",nstart," nend",nend
411 c read(inp,'(a)') pdbfile
412 c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
413 c open(ipdbin,file=pdbfile,status='old',err=33)
415 c 33 write (iout,'(a)') 'Error opening PDB file.'
418 c print *,'Begin reading pdb data'
420 c print *,'Finished reading pdb data'
421 c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
423 c itype_pdb(i)=itype(i)
426 c write (iout,'(a,i3)') 'nsup=',nsup
428 c if (nsup.le.(nct-nnt+1)) then
429 c do i=0,nct-nnt+1-nsup
430 c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
436 c & 'Error - sequences to be superposed do not match.'
439 c do i=0,nsup-(nct-nnt+1)
440 c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
442 c nstart_sup=nstart_sup+i
448 c & 'Error - sequences to be superposed do not match.'
451 c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
452 c & ' nstart_seq=',nstart_seq
456 write (iout,*) "molread: REFSTR",refstr
458 if (.not.pdbref) then
459 call read_angles(inp,*38)
461 38 write (iout,'(a)') 'Error reading reference structure.'
463 call mp_stopall(Error_Msg)
465 stop 'Error reading reference structure'
478 call contact(.true.,ncont_ref,icont_ref)
481 C write (iout,'(/a,i3,a)')
482 C & 'The chain contains',ns,' disulfide-bridging cysteines.'
483 write (iout,'(20i4)') (iss(i),i=1,ns)
485 write(iout,*)"Running with dynamic disulfide-bond formation"
487 write (iout,'(/a/)') 'Pre-formed links are:'
493 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
494 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
500 if (ns.gt.0.and.dyn_ss) then
504 forcon(i-nss)=forcon(i)
511 dyn_ss_mask(iss(i))=.true.
514 c Read distance restraints
515 if (constr_dist.gt.0) then
516 call read_dist_constr
521 c-----------------------------------------------------------------------------
522 logical function seq_comp(itypea,itypeb,length)
524 integer length,itypea(length),itypeb(length)
527 if (itypea(i).ne.itypeb(i)) then
535 c-----------------------------------------------------------------------------
536 subroutine read_bridge
537 C Read information about disulfide bridges.
540 include 'COMMON.IOUNITS'
543 include 'COMMON.INTERACT'
544 include 'COMMON.LOCAL'
545 include 'COMMON.NAMES'
546 include 'COMMON.CHAIN'
547 include 'COMMON.FFIELD'
548 include 'COMMON.SBRIDGE'
549 include 'COMMON.HEADER'
550 include 'COMMON.CONTROL'
551 include 'COMMON.TIME1'
553 include 'COMMON.INFO'
556 C Read bridging residues.
557 read (inp,*) ns,(iss(i),i=1,ns)
559 C Check whether the specified bridging residues are cystines.
561 if (itype(iss(i)).ne.1) then
562 write (iout,'(2a,i3,a)')
563 & 'Do you REALLY think that the residue ',
564 & restyp(itype(iss(i))),i,
565 & ' can form a disulfide bridge?!!!'
566 write (*,'(2a,i3,a)')
567 & 'Do you REALLY think that the residue ',
568 & restyp(itype(iss(i))),i,
569 & ' can form a disulfide bridge?!!!'
571 call mp_stopall(error_msg)
577 C Read preformed bridges.
579 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
582 C Check if the residues involved in bridges are in the specified list of
586 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
587 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
588 write (iout,'(a,i3,a)') 'Disulfide pair',i,
589 & ' contains residues present in other pairs.'
590 write (*,'(a,i3,a)') 'Disulfide pair',i,
591 & ' contains residues present in other pairs.'
593 call mp_stopall(error_msg)
600 if (ihpb(i).eq.iss(j)) goto 10
602 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
605 if (jhpb(i).eq.iss(j)) goto 20
607 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
620 c----------------------------------------------------------------------------
621 subroutine read_angles(kanal,*)
626 include 'COMMON.CHAIN'
627 include 'COMMON.IOUNITS'
629 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
630 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
631 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
632 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
634 theta(i)=deg2rad*theta(i)
635 phi(i)=deg2rad*phi(i)
636 alph(i)=deg2rad*alph(i)
637 omeg(i)=deg2rad*omeg(i)
642 c----------------------------------------------------------------------------
643 subroutine reada(rekord,lancuch,wartosc,default)
645 character*(*) rekord,lancuch
646 double precision wartosc,default
649 iread=index(rekord,lancuch)
654 iread=iread+ilen(lancuch)+1
655 read (rekord(iread:),*) wartosc
658 c----------------------------------------------------------------------------
659 subroutine multreada(rekord,lancuch,tablica,dim,default)
662 double precision tablica(dim),default
663 character*(*) rekord,lancuch
669 iread=index(rekord,lancuch)
670 if (iread.eq.0) return
671 iread=iread+ilen(lancuch)+1
672 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
675 c----------------------------------------------------------------------------
676 subroutine readi(rekord,lancuch,wartosc,default)
678 character*(*) rekord,lancuch
679 integer wartosc,default
682 iread=index(rekord,lancuch)
687 iread=iread+ilen(lancuch)+1
688 read (rekord(iread:),*) wartosc
691 C----------------------------------------------------------------------
692 subroutine multreadi(rekord,lancuch,tablica,dim,default)
695 integer tablica(dim),default
696 character*(*) rekord,lancuch
703 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
704 if (iread.eq.0) return
705 iread=iread+ilen(lancuch)+1
706 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
710 c----------------------------------------------------------------------------
711 subroutine card_concat(card)
713 include 'COMMON.IOUNITS'
715 character*80 karta,ucase
717 read (inp,'(a)') karta
720 do while (karta(80:80).eq.'&')
721 card=card(:ilen(card)+1)//karta(:79)
722 read (inp,'(a)') karta
725 card=card(:ilen(card)+1)//karta
728 c----------------------------------------------------------------------------
737 include 'COMMON.IOUNITS'
738 include 'COMMON.CONTROL'
739 integer lenpre,lenpot,ilen
741 character*16 cformat,cprint
743 integer lenint,lenout
744 call getenv('INPUT',prefix)
745 call getenv('OUTPUT',prefout)
746 call getenv('INTIN',prefintin)
747 call getenv('COORD',cformat)
748 call getenv('PRINTCOOR',cprint)
749 call getenv('SCRATCHDIR',scratchdir)
752 if (index(ucase(cformat),'CX').gt.0) then
759 lenint=ilen(prefintin)
760 C Get the names and open the input files
761 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
763 write (liczba,'(bz,i3.3)') me
764 outname=prefout(:lenout)//'_clust.out_'//liczba
766 outname=prefout(:lenout)//'_clust.out'
769 intinname=prefintin(:lenint)//'.bx'
770 else if (from_cx) then
771 intinname=prefintin(:lenint)//'.cx'
773 intinname=prefintin(:lenint)//'.int'
775 rmsname=prefintin(:lenint)//'.rms'
776 open (jplot,file=prefout(:ilen(prefout))//'.tex',
778 open (jrms,file=rmsname,status='unknown')
779 open(iout,file=outname,status='unknown')
780 C Get parameter filenames and open the parameter files.
781 call getenv('BONDPAR',bondname)
782 open (ibond,file=bondname,status='old')
783 call getenv('THETPAR',thetname)
784 open (ithep,file=thetname,status='old')
785 call getenv('ROTPAR',rotname)
786 open (irotam,file=rotname,status='old')
787 call getenv('TORPAR',torname)
788 open (itorp,file=torname,status='old')
789 call getenv('TORDPAR',tordname)
790 open (itordp,file=tordname,status='old')
791 call getenv('FOURIER',fouriername)
792 open (ifourier,file=fouriername,status='old')
793 call getenv('ELEPAR',elename)
794 open (ielep,file=elename,status='old')
795 call getenv('SIDEPAR',sidename)
796 open (isidep,file=sidename,status='old')
797 call getenv('SIDEP',sidepname)
798 open (isidep1,file=sidepname,status="old")
799 call getenv('SCCORPAR',sccorname)
800 open (isccor,file=sccorname,status="old")
801 call getenv('LIPTRANPAR',liptranname)
802 open (iliptranpar,file=liptranname,status='old')
805 C 8/9/01 In the newest version SCp interaction constants are read from a file
806 C Use -DOLDSCP to use hard-coded constants instead.
808 call getenv('SCPPAR',scpname)
809 open (iscpp,file=scpname,status='old')
813 subroutine read_dist_constr
814 implicit real*8 (a-h,o-z)
819 include 'COMMON.CONTROL'
820 include 'COMMON.CHAIN'
821 include 'COMMON.IOUNITS'
822 include 'COMMON.SBRIDGE'
823 integer ifrag_(2,100),ipair_(2,100)
824 double precision wfrag_(100),wpair_(100)
825 character*500 controlcard
826 logical lprn /.true./
827 write (iout,*) "Calling read_dist_constr"
828 C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
830 write(iout,*) "TU sie wywalam?"
831 call card_concat(controlcard)
832 write (iout,*) controlcard
834 call readi(controlcard,"NFRAG",nfrag_,0)
835 call readi(controlcard,"NPAIR",npair_,0)
836 call readi(controlcard,"NDIST",ndist_,0)
837 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
838 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
839 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
840 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
841 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
842 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
843 write (iout,*) "IFRAG"
845 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
847 write (iout,*) "IPAIR"
849 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
853 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
854 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
855 & ifrag_(2,i)=nstart_sup+nsup-1
856 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
858 if (wfrag_(i).gt.0.0d0) then
859 do j=ifrag_(1,i),ifrag_(2,i)-1
861 write (iout,*) "j",j," k",k
863 if (constr_dist.eq.1) then
868 forcon(nhpb)=wfrag_(i)
869 else if (constr_dist.eq.2) then
870 if (ddjk.le.dist_cut) then
875 forcon(nhpb)=wfrag_(i)
882 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
885 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
886 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
892 if (wpair_(i).gt.0.0d0) then
900 do j=ifrag_(1,ii),ifrag_(2,ii)
901 do k=ifrag_(1,jj),ifrag_(2,jj)
905 forcon(nhpb)=wpair_(i)
907 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
908 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
914 if (constr_dist.eq.11) then
915 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
916 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
917 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
918 C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
919 C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
921 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
923 if (forcon(nhpb+1).gt.0.0d0) then
925 if (ibecarb(i).gt.0) then
929 if (dhpb(nhpb).eq.0.0d0)
930 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
931 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
932 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
933 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)