3 C Read the parameters of the probability distributions of the virtual-bond
4 C valence angles and the side chains and energy parameters.
6 C Important! Energy-term weights ARE NOT read here; they are read from the
7 C main input file instead, because NO defaults have yet been set for these
10 implicit real*8 (a-h,o-z)
16 include 'COMMON.IOUNITS'
17 include 'COMMON.CHAIN'
18 include 'COMMON.INTERACT'
20 include 'COMMON.LOCAL'
21 include 'COMMON.TORSION'
22 include 'COMMON.SCCOR'
23 include 'COMMON.SCROT'
24 include 'COMMON.FFIELD'
25 include 'COMMON.NAMES'
26 include 'COMMON.SBRIDGE'
28 include 'COMMON.SETUP'
30 character*1 onelett(4) /"G","A","P","D"/
31 character*1 toronelet(-2:2) /"p","a","G","A","P"/
33 dimension blower(3,3,maxlob)
35 character*3 lancuch,ucase
37 C For printing parameters after they are read set the following in the UNRES
40 C setenv PRINT_PARM YES
42 C To print parameters in LaTeX format rather than as ASCII tables:
46 call getenv_loc("PRINT_PARM",lancuch)
47 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
48 call getenv_loc("LATEX",lancuch)
49 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
51 dwa16=2.0d0**(1.0d0/6.0d0)
53 C Assign virtual-bond length
58 c Read the virtual-bond parameters, masses, and moments of inertia
59 c and Stokes' radii of the peptide group and side chains
62 read (ibond,*) vbldp0,akp,mp,ip,pstok
65 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
70 dsc_inv(i)=1.0D0/dsc(i)
74 read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
76 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
77 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
82 dsc_inv(i)=1.0D0/dsc(i)
87 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
88 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
90 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
92 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
93 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
95 write (iout,'(13x,3f10.5)')
96 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
102 C Read the parameters of the probability distribution/energy expression
103 C of the virtual-bond valence angles theta
106 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
107 & (bthet(j,i,1,1),j=1,2)
108 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
109 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
110 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
114 athet(1,i,1,-1)=athet(1,i,1,1)
115 athet(2,i,1,-1)=athet(2,i,1,1)
116 bthet(1,i,1,-1)=-bthet(1,i,1,1)
117 bthet(2,i,1,-1)=-bthet(2,i,1,1)
118 athet(1,i,-1,1)=-athet(1,i,1,1)
119 athet(2,i,-1,1)=-athet(2,i,1,1)
120 bthet(1,i,-1,1)=bthet(1,i,1,1)
121 bthet(2,i,-1,1)=bthet(2,i,1,1)
125 athet(1,i,-1,-1)=athet(1,-i,1,1)
126 athet(2,i,-1,-1)=-athet(2,-i,1,1)
127 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
128 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
129 athet(1,i,-1,1)=athet(1,-i,1,1)
130 athet(2,i,-1,1)=-athet(2,-i,1,1)
131 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
132 bthet(2,i,-1,1)=bthet(2,-i,1,1)
133 athet(1,i,1,-1)=-athet(1,-i,1,1)
134 athet(2,i,1,-1)=athet(2,-i,1,1)
135 bthet(1,i,1,-1)=bthet(1,-i,1,1)
136 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
141 polthet(j,i)=polthet(j,-i)
144 gthet(j,i)=gthet(j,-i)
151 & 'Parameters of the virtual-bond valence angles:'
152 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
153 & ' ATHETA0 ',' A1 ',' A2 ',
156 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
157 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
159 write (iout,'(/a/9x,5a/79(1h-))')
160 & 'Parameters of the expression for sigma(theta_c):',
161 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
162 & ' ALPH3 ',' SIGMA0C '
164 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
165 & (polthet(j,i),j=0,3),sigc0(i)
167 write (iout,'(/a/9x,5a/79(1h-))')
168 & 'Parameters of the second gaussian:',
169 & ' THETA0 ',' SIGMA0 ',' G1 ',
172 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
173 & sig0(i),(gthet(j,i),j=1,3)
177 & 'Parameters of the virtual-bond valence angles:'
178 write (iout,'(/a/9x,5a/79(1h-))')
179 & 'Coefficients of expansion',
180 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
181 & ' b1*10^1 ',' b2*10^1 '
183 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
184 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
185 $ (10*bthet(j,i,1,1),j=1,2)
187 write (iout,'(/a/9x,5a/79(1h-))')
188 & 'Parameters of the expression for sigma(theta_c):',
189 & ' alpha0 ',' alph1 ',' alph2 ',
190 & ' alhp3 ',' sigma0c '
192 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
193 & (polthet(j,i),j=0,3),sigc0(i)
195 write (iout,'(/a/9x,5a/79(1h-))')
196 & 'Parameters of the second gaussian:',
197 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
200 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
201 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
207 C Read the parameters of Utheta determined from ab initio surfaces
208 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
210 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
211 & ntheterm3,nsingle,ndouble
212 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
213 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
219 aathet(l,i,j,k)=0.0d0
223 bbthet(m,l,i,j,k)=0.0d0
224 ccthet(m,l,i,j,k)=0.0d0
225 ddthet(m,l,i,j,k)=0.0d0
226 eethet(m,l,i,j,k)=0.0d0
232 ffthet(mm,m,l,i,j,k)=0.0d0
233 ggthet(mm,m,l,i,j,k)=0.0d0
243 read (ithep,'(3a)',end=111,err=111) res1,res2,res3
244 read (ithep,*,end=111,err=111) aa0thet(i,j,k)
245 read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
246 read (ithep,*,end=111,err=111)
247 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
248 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
249 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
250 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
251 read (ithep,*,end=111,err=111)
252 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
253 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
254 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
259 C For dummy ends assign glycine-type coefficients of theta-only terms; the
260 C coefficients of theta-and-gamma-dependent terms are zero.
265 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
266 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
268 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
269 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
272 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
274 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
277 C Control printout of the coefficients of virtual-bond-angle potentials
280 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
284 write (iout,'(//4a)')
285 & 'Type ',onelett(i),onelett(j),onelett(k)
286 write (iout,'(//a,10x,a)') " l","a[l]"
287 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
288 write (iout,'(i2,1pe15.5)')
289 & (l,aathet(l,i,j,k),l=1,ntheterm)
291 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
292 & "b",l,"c",l,"d",l,"e",l
294 write (iout,'(i2,4(1pe15.5))') m,
295 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
296 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
300 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
301 & "f+",l,"f-",l,"g+",l,"g-",l
304 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
305 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
306 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
315 write (2,*) "Start reading THETA_PDB"
317 read (ithep_pdb,*,err=111,end=111) a0thet(i),
318 & (athet(j,i,1,1),j=1,2),
319 & (bthet(j,i,1,1),j=1,2)
320 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
321 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
322 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
324 write (2,*) "End reading THETA_PDB"
331 C Read the parameters of the probability distribution/energy expression
332 C of the side chains.
335 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
339 dsc_inv(i)=1.0D0/dsc(i)
350 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
351 & ((blower(k,l,1),l=1,k),k=1,3)
352 censc(1,1,-i)=censc(1,1,i)
353 censc(2,1,-i)=censc(2,1,i)
354 censc(3,1,-i)=-censc(3,1,i)
357 read (irotam,*,end=112,err=112) bsc(j,i)
358 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
359 & ((blower(k,l,j),l=1,k),k=1,3)
360 censc(1,j,-i)=censc(1,j,i)
361 censc(2,j,-i)=censc(2,j,i)
362 censc(3,j,-i)=-censc(3,j,i)
363 C BSC is amplitude of Gaussian
370 akl=akl+blower(k,m,j)*blower(l,m,j)
374 if (((k.eq.3).and.(l.ne.3))
375 & .or.((l.eq.3).and.(k.ne.3))) then
376 gaussc(k,l,j,-i)=-akl
377 gaussc(l,k,j,-i)=-akl
389 write (iout,'(/a)') 'Parameters of side-chain local geometry'
394 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
395 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
396 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
397 & 'log h',(bsc(j,i),j=1,nlobi)
398 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
399 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
401 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
402 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
405 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
406 write (iout,'(a,f10.4,4(16x,f10.4))')
407 & 'Center ',(bsc(j,i),j=1,nlobi)
408 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
417 C Read scrot parameters for potentials determined from all-atom AM1 calculations
418 C added by Urszula Kozlowska 07/11/2007
421 read (irotam,*,end=112,err=112)
423 read (irotam,*,end=112,err=112)
426 read(irotam,*,end=112,err=112) sc_parmin(j,i)
431 C Read the parameters of the probability distribution/energy expression
432 C of the side chains.
435 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
439 dsc_inv(i)=1.0D0/dsc(i)
450 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
451 & ((blower(k,l,1),l=1,k),k=1,3)
453 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
454 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
455 & ((blower(k,l,j),l=1,k),k=1,3)
462 akl=akl+blower(k,m,j)*blower(l,m,j)
477 C Read torsional parameters in old format
479 read (itorp,*,end=113,err=113) ntortyp,nterm_old
480 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
481 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
486 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
492 write (iout,'(/a/)') 'Torsional constants:'
495 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
496 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
502 C Read torsional parameters
504 read (itorp,*,end=113,err=113) ntortyp
505 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
508 itortyp(i)=-itortyp(-i)
510 c write (iout,*) 'ntortyp',ntortyp
512 do j=-ntortyp+1,ntortyp-1
513 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
515 nterm(-i,-j,iblock)=nterm(i,j,iblock)
516 nlor(-i,-j,iblock)=nlor(i,j,iblock)
519 c write (iout,*) i,j,nterm(i,j,iblock),nlor(i,j,iblock)
520 do k=1,nterm(i,j,iblock)
521 c write(iout,*) 'wchodze',k
522 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
524 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
525 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
526 v0ij=v0ij+si*v1(k,i,j,iblock)
528 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
529 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
530 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
532 do k=1,nlor(i,j,iblock)
533 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
534 & vlor2(k,i,j),vlor3(k,i,j)
535 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
538 v0(-i,-j,iblock)=v0ij
544 write (iout,'(/a/)') 'Torsional constants:'
547 write (iout,*) 'ityp',i,' jtyp',j
548 write (iout,*) 'Fourier constants'
549 do k=1,nterm(i,j,iblock)
550 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
553 write (iout,*) 'Lorenz constants'
554 do k=1,nlor(i,j,iblock)
555 write (iout,'(3(1pe15.5))')
556 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
562 C 6/23/01 Read parameters for double torsionals
566 do j=-ntortyp+1,ntortyp-1
567 do k=-ntortyp+1,ntortyp-1
568 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
569 write (iout,*) "OK onelett",
572 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
573 & .or. t3.ne.toronelet(k)) then
574 write (iout,*) "Error in double torsional parameter file",
577 call MPI_Finalize(Ierror)
579 stop "Error in double torsional parameter file"
581 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
582 & ntermd_2(i,j,k,iblock)
583 write (iout,*) ntermd_1(i,j,k,iblock),
584 & ntermd_2(i,j,k,iblock)
585 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
586 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
587 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
588 & ntermd_1(i,j,k,iblock))
589 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
590 & ntermd_1(i,j,k,iblock))
591 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
592 & ntermd_1(i,j,k,iblock))
593 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
594 & ntermd_1(i,j,k,iblock))
595 C Matrix of D parameters for one dimesional foureir series
596 do l=1,ntermd_1(i,j,k,iblock)
597 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
598 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
599 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
600 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
601 c write(iout,*) "whcodze" ,
602 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
604 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
605 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
606 & v2s(m,l,i,j,k,iblock),
607 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
608 C Matrix of D parameters for two dimesional fourier series
609 do l=1,ntermd_2(i,j,k,iblock)
611 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
612 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
613 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
614 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
623 write (iout,*) 'Constants for double torsionals'
626 do j=-ntortyp,ntortyp
627 do k=-ntortyp,ntortyp
628 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
629 & ' nsingle',ntermd_1(i,j,k,iblock),
630 & ' ndouble',ntermd_2(i,j,k,iblock)
632 write (iout,*) 'Single angles:'
633 do l=1,ntermd_1(i,j,k,iblock)
634 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
635 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
636 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
639 write (iout,*) 'Pairs of angles:'
640 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
641 do l=1,ntermd_2(i,j,k,iblock)
642 write (iout,'(i5,20f10.5)')
643 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
646 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
647 do l=1,ntermd_2(i,j,k,iblock)
648 write (iout,'(i5,20f10.5)')
649 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
650 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
660 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
661 C correlation energies.
663 read (isccor,*,end=119,err=119) nterm_sccor
668 read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
675 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
678 write (iout,*) 'ityp',i,' jtyp',j
680 write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
686 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
687 C interaction energy of the Gly, Ala, and Pro prototypes.
691 write (iout,*) "Coefficients of the cumulants"
693 read (ifourier,*) nloctyp
695 read (ifourier,*,end=115,err=115)
696 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
698 write (iout,*) 'Type',i
699 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
740 c Ctilde(1,1,i)=0.0d0
741 c Ctilde(1,2,i)=0.0d0
742 c Ctilde(2,1,i)=0.0d0
743 c Ctilde(2,2,i)=0.0d0
764 c Dtilde(1,1,i)=0.0d0
765 c Dtilde(1,2,i)=0.0d0
766 c Dtilde(2,1,i)=0.0d0
767 c Dtilde(2,2,i)=0.0d0
768 EE(1,1,i)= b(10)+b(11)
769 EE(2,2,i)=-b(10)+b(11)
770 EE(2,1,i)= b(12)-b(13)
771 EE(1,2,i)= b(12)+b(13)
772 EE(1,1,-i)= b(10)+b(11)
773 EE(2,2,-i)=-b(10)+b(11)
774 EE(2,1,-i)=-b(12)+b(13)
775 EE(1,2,-i)=-b(12)-b(13)
780 c ee(2,1,i)=ee(1,2,i)
784 write (iout,*) 'Type',i
786 write(iout,*) B1(1,i),B1(2,i)
788 write(iout,*) B2(1,i),B2(2,i)
791 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
795 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
799 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
804 C Read electrostatic-interaction parameters
808 write (iout,'(/a)') 'Electrostatic interaction constants:'
809 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
810 & 'IT','JT','APP','BPP','AEL6','AEL3'
812 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
813 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
814 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
815 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
820 app (i,j)=epp(i,j)*rri*rri
821 bpp (i,j)=-2.0D0*epp(i,j)*rri
822 ael6(i,j)=elpp6(i,j)*4.2D0**6
823 ael3(i,j)=elpp3(i,j)*4.2D0**3
824 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
825 & ael6(i,j),ael3(i,j)
829 C Read side-chain interaction parameters.
831 read (isidep,*,end=117,err=117) ipot,expon
832 if (ipot.lt.1 .or. ipot.gt.5) then
833 write (iout,'(2a)') 'Error while reading SC interaction',
834 & 'potential file - unknown potential type.'
836 call MPI_Finalize(Ierror)
842 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
843 & ', exponents are ',expon,2*expon
844 goto (10,20,30,30,40) ipot
845 C----------------------- LJ potential ---------------------------------
846 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
847 & (sigma0(i),i=1,ntyp)
849 write (iout,'(/a/)') 'Parameters of the LJ potential:'
850 write (iout,'(a/)') 'The epsilon array:'
851 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
852 write (iout,'(/a)') 'One-body parameters:'
853 write (iout,'(a,4x,a)') 'residue','sigma'
854 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
857 C----------------------- LJK potential --------------------------------
858 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
859 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
861 write (iout,'(/a/)') 'Parameters of the LJK potential:'
862 write (iout,'(a/)') 'The epsilon array:'
863 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
864 write (iout,'(/a)') 'One-body parameters:'
865 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
866 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
870 C---------------------- GB or BP potential -----------------------------
871 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
872 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
874 C For the GB potential convert sigma'**2 into chi'
877 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
881 write (iout,'(/a/)') 'Parameters of the BP potential:'
882 write (iout,'(a/)') 'The epsilon array:'
883 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
884 write (iout,'(/a)') 'One-body parameters:'
885 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
887 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
888 & chip(i),alp(i),i=1,ntyp)
891 C--------------------- GBV potential -----------------------------------
892 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
893 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
894 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
896 write (iout,'(/a/)') 'Parameters of the GBV potential:'
897 write (iout,'(a/)') 'The epsilon array:'
898 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
899 write (iout,'(/a)') 'One-body parameters:'
900 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
901 & 's||/s_|_^2',' chip ',' alph '
902 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
903 & sigii(i),chip(i),alp(i),i=1,ntyp)
907 C-----------------------------------------------------------------------
908 C Calculate the "working" parameters of SC interactions.
916 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
917 sigma(j,i)=sigma(i,j)
918 rs0(i,j)=dwa16*sigma(i,j)
922 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
923 & 'Working parameters of the SC interactions:',
924 & ' a ',' b ',' augm ',' sigma ',' r0 ',
929 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
938 sigeps=dsign(1.0D0,epsij)
940 aa(i,j)=epsij*rrij*rrij
941 bb(i,j)=-sigeps*epsij*rrij
949 ratsig1=sigt2sq/sigt1sq
950 ratsig2=1.0D0/ratsig1
951 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
952 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
953 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
957 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
958 sigmaii(i,j)=rsum_max
959 sigmaii(j,i)=rsum_max
961 c sigmaii(i,j)=r0(i,j)
962 c sigmaii(j,i)=r0(i,j)
964 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
965 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
966 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
967 augm(i,j)=epsij*r_augm**(2*expon)
968 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
975 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
976 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
977 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
983 C Define the SC-p interaction constants (hard-coded; old style)
986 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
988 c aad(i,1)=0.3D0*4.0D0**12
989 C Following line for constants currently implemented
990 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
991 aad(i,1)=1.5D0*4.0D0**12
992 c aad(i,1)=0.17D0*5.6D0**12
994 C "Soft" SC-p repulsion
996 C Following line for constants currently implemented
997 c aad(i,1)=0.3D0*4.0D0**6
998 C "Hard" SC-p repulsion
999 bad(i,1)=3.0D0*4.0D0**6
1000 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1009 C 8/9/01 Read the SC-p interaction constants from file
1012 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1015 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1016 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1017 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1018 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1022 write (iout,*) "Parameters of SC-p interactions:"
1024 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1025 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1030 C Define the constants of the disulfide bridge
1034 c Old arbitrary potential - commented out.
1039 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1040 c energy surface of diethyl disulfide.
1041 c A. Liwo and U. Kozlowska, 11/24/03
1058 write (iout,'(/a)') "Disulfide bridge parameters:"
1059 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1060 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1061 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1062 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1066 111 write (iout,*) "Error reading bending energy parameters."
1068 112 write (iout,*) "Error reading rotamer energy parameters."
1070 113 write (iout,*) "Error reading torsional energy parameters."
1072 114 write (iout,*) "Error reading double torsional energy parameters."
1075 & "Error reading cumulant (multibody energy) parameters."
1077 116 write (iout,*) "Error reading electrostatic energy parameters."
1079 117 write (iout,*) "Error reading side chain interaction parameters."
1081 118 write (iout,*) "Error reading SCp interaction parameters."
1083 119 write (iout,*) "Error reading SCCOR parameters"
1086 call MPI_Finalize(Ierror)
1093 subroutine getenv_loc(var, val)
1094 character(*) var, val
1097 character(2000) line
1100 open (196,file='env',status='old',readonly,shared)
1102 c write(*,*)'looking for ',var
1103 10 read(196,*,err=11,end=11)line
1104 iread=index(line,var)
1105 c write(*,*)iread,' ',var,' ',line
1106 if (iread.eq.0) go to 10
1107 c write(*,*)'---> ',line
1113 iread=iread+ilen(var)+1
1114 read (line(iread:),*,err=12,end=12) val
1115 c write(*,*)'OK: ',var,' = ',val
1121 #elif (defined CRAY)
1122 integer lennam,lenval,ierror
1124 c getenv using a POSIX call, useful on the T3D
1125 c Sept 1996, comment out error check on advice of H. Pritchard
1128 if(lennam.le.0) stop '--error calling getenv--'
1129 call pxfgetenv(var,lennam,val,lenval,ierror)
1130 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1132 call getenv(var,val)