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),(athet(j,i),j=1,2),
319 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
320 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
321 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
324 write (2,*) "End reading THETA_PDB"
330 C Read the parameters of the probability distribution/energy expression
331 C of the side chains.
334 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
338 dsc_inv(i)=1.0D0/dsc(i)
349 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
350 & ((blower(k,l,1),l=1,k),k=1,3)
351 censc(1,1,-i)=censc(1,1,i)
352 censc(2,1,-i)=censc(2,1,i)
353 censc(3,1,-i)=-censc(3,1,i)
356 read (irotam,*,end=112,err=112) bsc(j,i)
357 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
358 & ((blower(k,l,j),l=1,k),k=1,3)
359 censc(1,j,-i)=censc(1,j,i)
360 censc(2,j,-i)=censc(2,j,i)
361 censc(3,j,-i)=-censc(3,j,i)
362 C BSC is amplitude of Gaussian
369 akl=akl+blower(k,m,j)*blower(l,m,j)
373 if (((k.eq.3).and.(l.ne.3))
374 & .or.((l.eq.3).and.(k.ne.3))) then
375 gaussc(k,l,j,-i)=-akl
376 gaussc(l,k,j,-i)=-akl
388 write (iout,'(/a)') 'Parameters of side-chain local geometry'
393 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
394 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
395 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
396 & 'log h',(bsc(j,i),j=1,nlobi)
397 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
398 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
400 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
401 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
404 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
405 write (iout,'(a,f10.4,4(16x,f10.4))')
406 & 'Center ',(bsc(j,i),j=1,nlobi)
407 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
416 C Read scrot parameters for potentials determined from all-atom AM1 calculations
417 C added by Urszula Kozlowska 07/11/2007
420 read (irotam,*,end=112,err=112)
422 read (irotam,*,end=112,err=112)
425 read(irotam,*,end=112,err=112) sc_parmin(j,i)
430 C Read the parameters of the probability distribution/energy expression
431 C of the side chains.
434 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
438 dsc_inv(i)=1.0D0/dsc(i)
449 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
450 & ((blower(k,l,1),l=1,k),k=1,3)
452 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
453 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
454 & ((blower(k,l,j),l=1,k),k=1,3)
461 akl=akl+blower(k,m,j)*blower(l,m,j)
476 C Read torsional parameters in old format
478 read (itorp,*,end=113,err=113) ntortyp,nterm_old
479 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
480 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
485 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
491 write (iout,'(/a/)') 'Torsional constants:'
494 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
495 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
501 C Read torsional parameters
503 read (itorp,*,end=113,err=113) ntortyp
504 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
507 itortyp(i)=-itortyp(-i)
509 c write (iout,*) 'ntortyp',ntortyp
511 do j=-ntortyp+1,ntortyp-1
512 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
514 nterm(-i,-j,iblock)=nterm(i,j,iblock)
515 nlor(-i,-j,iblock)=nlor(i,j,iblock)
518 do k=1,nterm(i,j,iblock)
519 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
521 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
522 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
523 v0ij=v0ij+si*v1(k,i,j,iblock)
525 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
526 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
527 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
529 do k=1,nlor(i,j,iblock)
530 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
531 & vlor2(k,i,j),vlor3(k,i,j)
532 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
535 v0(-i,-j,iblock)=v0ij
541 write (iout,'(/a/)') 'Torsional constants:'
544 write (iout,*) 'ityp',i,' jtyp',j
545 write (iout,*) 'Fourier constants'
546 do k=1,nterm(i,j,iblock)
547 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
550 write (iout,*) 'Lorenz constants'
551 do k=1,nlor(i,j,iblock)
552 write (iout,'(3(1pe15.5))')
553 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
559 C 6/23/01 Read parameters for double torsionals
563 do j=-ntortyp+1,ntortyp-1
564 do k=-ntortyp+1,ntortyp-1
565 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
566 c write (iout,*) "OK onelett",
569 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
570 & .or. t3.ne.toronelet(k)) then
571 write (iout,*) "Error in double torsional parameter file",
574 call MPI_Finalize(Ierror)
576 stop "Error in double torsional parameter file"
578 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
579 & ntermd_2(i,j,k,iblock)
580 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
581 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
582 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
583 & ntermd_1(i,j,k,iblock))
584 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
585 & ntermd_1(i,j,k,iblock))
586 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
587 & ntermd_1(i,j,k,iblock))
588 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
589 & ntermd_1(i,j,k,iblock))
590 C Martix of D parameters for one dimesional foureir series
591 do l=1,ntermd_1(i,j,k,iblock)
592 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
593 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
594 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
595 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
596 c write(iout,*) "whcodze" ,
597 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
599 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
600 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
601 & v2s(m,l,i,j,k,iblock),
602 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
603 C Martix of D parameters for two dimesional fourier series
604 do l=1,ntermd_2(i,j,k,iblock)
606 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
607 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
608 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
609 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
618 write (iout,*) 'Constants for double torsionals'
621 do j=-ntortyp,ntortyp
622 do k=-ntortyp,ntortyp
623 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
624 & ' nsingle',ntermd_1(i,j,k,iblock),
625 & ' ndouble',ntermd_2(i,j,k,iblock)
627 write (iout,*) 'Single angles:'
628 do l=1,ntermd_1(i,j,k,iblock)
629 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
630 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
631 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
632 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
635 write (iout,*) 'Pairs of angles:'
636 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
637 do l=1,ntermd_2(i,j,k,iblock)
638 write (iout,'(i5,20f10.5)')
639 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
642 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
643 do l=1,ntermd_2(i,j,k,iblock)
644 write (iout,'(i5,20f10.5)')
645 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
646 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
655 C Read of Side-chain backbone correlation parameters
656 C Modified 11 May 2012 by Adasko
659 read (isccor,*,end=113,err=113) nsccortyp
661 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
663 isccortyp(i)=-isccortyp(-i)
665 iscprol=isccortyp(20)
666 c write (iout,*) 'ntortyp',ntortyp
668 cc maxinter is maximum interaction sites
672 read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
676 do k=1,nterm_sccor(i,j)
677 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
679 if (j.eq.iscprol) then
680 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
681 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
682 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
683 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
684 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
685 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
686 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
687 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
689 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
690 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
691 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
692 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
693 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
694 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
696 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
699 do k=1,nlor_sccor(i,j)
700 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
701 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
702 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
703 &(1+vlor3sccor(k,i,j)**2)
705 v0sccor(i,j)=v0ijsccor
711 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
712 c write (iout,*) 'ntortyp',ntortyp
714 cc maxinter is maximum interaction sites
718 read (isccor,*,end=113,err=113)
719 & nterm_sccor(i,j),nlor_sccor(i,j)
723 do k=1,nterm_sccor(i,j)
724 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
726 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
729 do k=1,nlor_sccor(i,j)
730 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
731 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
732 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
733 &(1+vlor3sccor(k,i,j)**2)
735 v0sccor(i,j)=v0ijsccor
743 write (iout,'(/a/)') 'Torsional constants:'
746 write (iout,*) 'ityp',i,' jtyp',j
747 write (iout,*) 'Fourier constants'
748 do k=1,nterm_sccor(i,j)
749 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
751 write (iout,*) 'Lorenz constants'
752 do k=1,nlor_sccor(i,j)
753 write (iout,'(3(1pe15.5))')
754 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
761 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
762 C interaction energy of the Gly, Ala, and Pro prototypes.
766 write (iout,*) "Coefficients of the cumulants"
768 read (ifourier,*) nloctyp
770 read (ifourier,*,end=115,err=115)
771 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
773 write (iout,*) 'Type',i
774 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
816 c Ctilde(1,1,i)=0.0d0
817 c Ctilde(1,2,i)=0.0d0
818 c Ctilde(2,1,i)=0.0d0
819 c Ctilde(2,2,i)=0.0d0
841 c Dtilde(1,1,i)=0.0d0
842 c Dtilde(1,2,i)=0.0d0
843 c Dtilde(2,1,i)=0.0d0
844 c Dtilde(2,2,i)=0.0d0
845 EE(1,1,i)= b(10)+b(11)
846 EE(2,2,i)=-b(10)+b(11)
847 EE(2,1,i)= b(12)-b(13)
848 EE(1,2,i)= b(12)+b(13)
849 EE(1,1,-i)= b(10)+b(11)
850 EE(2,2,-i)=-b(10)+b(11)
851 EE(2,1,-i)=-b(12)+b(13)
852 EE(1,2,-i)=-b(12)-b(13)
858 c ee(2,1,i)=ee(1,2,i)
862 write (iout,*) 'Type',i
864 write(iout,*) B1(1,i),B1(2,i)
866 write(iout,*) B2(1,i),B2(2,i)
869 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
873 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
877 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
882 C Read electrostatic-interaction parameters
886 write (iout,'(/a)') 'Electrostatic interaction constants:'
887 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
888 & 'IT','JT','APP','BPP','AEL6','AEL3'
890 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
891 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
892 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
893 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
898 app (i,j)=epp(i,j)*rri*rri
899 bpp (i,j)=-2.0D0*epp(i,j)*rri
900 ael6(i,j)=elpp6(i,j)*4.2D0**6
901 ael3(i,j)=elpp3(i,j)*4.2D0**3
902 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
903 & ael6(i,j),ael3(i,j)
907 C Read side-chain interaction parameters.
909 read (isidep,*,end=117,err=117) ipot,expon
910 if (ipot.lt.1 .or. ipot.gt.5) then
911 write (iout,'(2a)') 'Error while reading SC interaction',
912 & 'potential file - unknown potential type.'
914 call MPI_Finalize(Ierror)
920 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
921 & ', exponents are ',expon,2*expon
922 goto (10,20,30,30,40) ipot
923 C----------------------- LJ potential ---------------------------------
924 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
925 & (sigma0(i),i=1,ntyp)
927 write (iout,'(/a/)') 'Parameters of the LJ potential:'
928 write (iout,'(a/)') 'The epsilon array:'
929 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
930 write (iout,'(/a)') 'One-body parameters:'
931 write (iout,'(a,4x,a)') 'residue','sigma'
932 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
935 C----------------------- LJK potential --------------------------------
936 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
937 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
939 write (iout,'(/a/)') 'Parameters of the LJK potential:'
940 write (iout,'(a/)') 'The epsilon array:'
941 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
942 write (iout,'(/a)') 'One-body parameters:'
943 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
944 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
948 C---------------------- GB or BP potential -----------------------------
949 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
950 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
952 C For the GB potential convert sigma'**2 into chi'
955 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
959 write (iout,'(/a/)') 'Parameters of the BP potential:'
960 write (iout,'(a/)') 'The epsilon array:'
961 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
962 write (iout,'(/a)') 'One-body parameters:'
963 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
965 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
966 & chip(i),alp(i),i=1,ntyp)
969 C--------------------- GBV potential -----------------------------------
970 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
971 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
972 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
974 write (iout,'(/a/)') 'Parameters of the GBV potential:'
975 write (iout,'(a/)') 'The epsilon array:'
976 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
977 write (iout,'(/a)') 'One-body parameters:'
978 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
979 & 's||/s_|_^2',' chip ',' alph '
980 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
981 & sigii(i),chip(i),alp(i),i=1,ntyp)
985 C-----------------------------------------------------------------------
986 C Calculate the "working" parameters of SC interactions.
994 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
995 sigma(j,i)=sigma(i,j)
996 rs0(i,j)=dwa16*sigma(i,j)
1000 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1001 & 'Working parameters of the SC interactions:',
1002 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1007 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1016 sigeps=dsign(1.0D0,epsij)
1018 aa(i,j)=epsij*rrij*rrij
1019 bb(i,j)=-sigeps*epsij*rrij
1023 sigt1sq=sigma0(i)**2
1024 sigt2sq=sigma0(j)**2
1027 ratsig1=sigt2sq/sigt1sq
1028 ratsig2=1.0D0/ratsig1
1029 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1030 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1031 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1035 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1036 sigmaii(i,j)=rsum_max
1037 sigmaii(j,i)=rsum_max
1039 c sigmaii(i,j)=r0(i,j)
1040 c sigmaii(j,i)=r0(i,j)
1042 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1043 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1044 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1045 augm(i,j)=epsij*r_augm**(2*expon)
1046 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1053 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1054 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1055 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1061 C Define the SC-p interaction constants (hard-coded; old style)
1064 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1066 c aad(i,1)=0.3D0*4.0D0**12
1067 C Following line for constants currently implemented
1068 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1069 aad(i,1)=1.5D0*4.0D0**12
1070 c aad(i,1)=0.17D0*5.6D0**12
1072 C "Soft" SC-p repulsion
1074 C Following line for constants currently implemented
1075 c aad(i,1)=0.3D0*4.0D0**6
1076 C "Hard" SC-p repulsion
1077 bad(i,1)=3.0D0*4.0D0**6
1078 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1087 C 8/9/01 Read the SC-p interaction constants from file
1090 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1093 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1094 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1095 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1096 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1100 write (iout,*) "Parameters of SC-p interactions:"
1102 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1103 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1108 C Define the constants of the disulfide bridge
1112 c Old arbitrary potential - commented out.
1117 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1118 c energy surface of diethyl disulfide.
1119 c A. Liwo and U. Kozlowska, 11/24/03
1136 write (iout,'(/a)') "Disulfide bridge parameters:"
1137 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1138 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1139 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1140 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1144 111 write (iout,*) "Error reading bending energy parameters."
1146 112 write (iout,*) "Error reading rotamer energy parameters."
1148 113 write (iout,*) "Error reading torsional energy parameters."
1150 114 write (iout,*) "Error reading double torsional energy parameters."
1153 & "Error reading cumulant (multibody energy) parameters."
1155 116 write (iout,*) "Error reading electrostatic energy parameters."
1157 117 write (iout,*) "Error reading side chain interaction parameters."
1159 118 write (iout,*) "Error reading SCp interaction parameters."
1161 119 write (iout,*) "Error reading SCCOR parameters"
1164 call MPI_Finalize(Ierror)
1171 subroutine getenv_loc(var, val)
1172 character(*) var, val
1175 character(2000) line
1178 open (196,file='env',status='old',readonly,shared)
1180 c write(*,*)'looking for ',var
1181 10 read(196,*,err=11,end=11)line
1182 iread=index(line,var)
1183 c write(*,*)iread,' ',var,' ',line
1184 if (iread.eq.0) go to 10
1185 c write(*,*)'---> ',line
1191 iread=iread+ilen(var)+1
1192 read (line(iread:),*,err=12,end=12) val
1193 c write(*,*)'OK: ',var,' = ',val
1199 #elif (defined CRAY)
1200 integer lennam,lenval,ierror
1202 c getenv using a POSIX call, useful on the T3D
1203 c Sept 1996, comment out error check on advice of H. Pritchard
1206 if(lennam.le.0) stop '--error calling getenv--'
1207 call pxfgetenv(var,lennam,val,lenval,ierror)
1208 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1210 call getenv(var,val)