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)
323 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,ntortyp
512 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
514 nterm(-i,-j,iblock)=nterm(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'
547 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
549 write (iout,*) 'Lorenz constants'
551 write (iout,'(3(1pe15.5))')
552 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
558 C 6/23/01 Read parameters for double torsionals
562 do j=-ntortyp+1,ntortyp-1
563 do k=-ntortyp+1,ntortyp-1
564 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
565 write (iout,*) "OK onelett",
568 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
569 & .or. t3.ne.toronelet(k)) then
570 write (iout,*) "Error in double torsional parameter file",
573 call MPI_Finalize(Ierror)
575 stop "Error in double torsional parameter file"
577 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
578 & ntermd_2(i,j,k,iblock)
579 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
580 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
581 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
582 & ntermd_1(i,j,k,iblock))
583 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
584 & ntermd_1(i,j,k,iblock))
585 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
586 & ntermd_1(i,j,k,iblock))
587 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
588 & ntermd_1(i,j,k,iblock))
589 C Matrix of D parameters for one dimesional foureir series
590 do l=1,ntermd_1(i,j,k,iblock)
591 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
592 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
593 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
594 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
595 c write(iout,*) "whcodze" ,
596 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
598 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
599 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
600 & v2s(m,l,i,j,k,iblock),
601 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
602 C Matrix of D parameters for two dimesional fourier series
603 do l=1,ntermd_2(i,j,k,iblock)
605 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
606 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
607 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
608 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
617 write (iout,*) 'Constants for double torsionals'
620 do j=-ntortyp,ntortyp
621 do k=-ntortyp,ntortyp
622 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
623 & ' nsingle',ntermd_1(i,j,k,iblock),
624 & ' ndouble',ntermd_2(i,j,k,iblock)
626 write (iout,*) 'Single angles:'
627 do l=1,ntermd_1(i,j,k,iblock)
628 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
629 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
630 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
633 write (iout,*) 'Pairs of angles:'
634 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
635 do l=1,ntermd_2(i,j,k,iblock)
636 write (iout,'(i5,20f10.5)')
637 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
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,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
644 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
654 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
655 C correlation energies.
657 read (isccor,*,end=119,err=119) nterm_sccor
662 read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
669 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
672 write (iout,*) 'ityp',i,' jtyp',j
674 write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
680 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
681 C interaction energy of the Gly, Ala, and Pro prototypes.
685 write (iout,*) "Coefficients of the cumulants"
687 read (ifourier,*) nloctyp
690 read (ifourier,*,end=115,err=115)
691 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
693 write (iout,*) 'Type',i
694 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
735 c Ctilde(1,1,i)=0.0d0
736 c Ctilde(1,2,i)=0.0d0
737 c Ctilde(2,1,i)=0.0d0
738 c Ctilde(2,2,i)=0.0d0
759 c Dtilde(1,1,i)=0.0d0
760 c Dtilde(1,2,i)=0.0d0
761 c Dtilde(2,1,i)=0.0d0
762 c Dtilde(2,2,i)=0.0d0
763 EE(1,1,i)= b(10)+b(11)
764 EE(2,2,i)=-b(10)+b(11)
765 EE(2,1,i)= b(12)-b(13)
766 EE(1,2,i)= b(12)+b(13)
767 EE(1,1,-i)= b(10)+b(11)
768 EE(2,2,-i)=-b(10)+b(11)
769 EE(2,1,-i)=-b(12)+b(13)
770 EE(1,2,-i)=-b(12)-b(13)
775 c ee(2,1,i)=ee(1,2,i)
779 write (iout,*) 'Type',i
781 write(iout,*) B1(1,i),B1(2,i)
783 write(iout,*) B2(1,i),B2(2,i)
786 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
790 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
794 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
799 C Read electrostatic-interaction parameters
803 write (iout,'(/a)') 'Electrostatic interaction constants:'
804 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
805 & 'IT','JT','APP','BPP','AEL6','AEL3'
807 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
808 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
809 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
810 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
815 app (i,j)=epp(i,j)*rri*rri
816 bpp (i,j)=-2.0D0*epp(i,j)*rri
817 ael6(i,j)=elpp6(i,j)*4.2D0**6
818 ael3(i,j)=elpp3(i,j)*4.2D0**3
819 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
820 & ael6(i,j),ael3(i,j)
824 C Read side-chain interaction parameters.
826 read (isidep,*,end=117,err=117) ipot,expon
827 if (ipot.lt.1 .or. ipot.gt.5) then
828 write (iout,'(2a)') 'Error while reading SC interaction',
829 & 'potential file - unknown potential type.'
831 call MPI_Finalize(Ierror)
837 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
838 & ', exponents are ',expon,2*expon
839 goto (10,20,30,30,40) ipot
840 C----------------------- LJ potential ---------------------------------
841 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
842 & (sigma0(i),i=1,ntyp)
844 write (iout,'(/a/)') 'Parameters of the LJ potential:'
845 write (iout,'(a/)') 'The epsilon array:'
846 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
847 write (iout,'(/a)') 'One-body parameters:'
848 write (iout,'(a,4x,a)') 'residue','sigma'
849 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
852 C----------------------- LJK potential --------------------------------
853 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
854 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
856 write (iout,'(/a/)') 'Parameters of the LJK potential:'
857 write (iout,'(a/)') 'The epsilon array:'
858 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
859 write (iout,'(/a)') 'One-body parameters:'
860 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
861 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
865 C---------------------- GB or BP potential -----------------------------
866 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
867 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
869 C For the GB potential convert sigma'**2 into chi'
872 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
876 write (iout,'(/a/)') 'Parameters of the BP potential:'
877 write (iout,'(a/)') 'The epsilon array:'
878 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
879 write (iout,'(/a)') 'One-body parameters:'
880 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
882 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
883 & chip(i),alp(i),i=1,ntyp)
886 C--------------------- GBV potential -----------------------------------
887 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
888 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
889 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
891 write (iout,'(/a/)') 'Parameters of the GBV potential:'
892 write (iout,'(a/)') 'The epsilon array:'
893 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
894 write (iout,'(/a)') 'One-body parameters:'
895 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
896 & 's||/s_|_^2',' chip ',' alph '
897 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
898 & sigii(i),chip(i),alp(i),i=1,ntyp)
902 C-----------------------------------------------------------------------
903 C Calculate the "working" parameters of SC interactions.
911 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
912 sigma(j,i)=sigma(i,j)
913 rs0(i,j)=dwa16*sigma(i,j)
917 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
918 & 'Working parameters of the SC interactions:',
919 & ' a ',' b ',' augm ',' sigma ',' r0 ',
924 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
933 sigeps=dsign(1.0D0,epsij)
935 aa(i,j)=epsij*rrij*rrij
936 bb(i,j)=-sigeps*epsij*rrij
944 ratsig1=sigt2sq/sigt1sq
945 ratsig2=1.0D0/ratsig1
946 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
947 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
948 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
952 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
953 sigmaii(i,j)=rsum_max
954 sigmaii(j,i)=rsum_max
956 c sigmaii(i,j)=r0(i,j)
957 c sigmaii(j,i)=r0(i,j)
959 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
960 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
961 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
962 augm(i,j)=epsij*r_augm**(2*expon)
963 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
970 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
971 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
972 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
978 C Define the SC-p interaction constants (hard-coded; old style)
981 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
983 c aad(i,1)=0.3D0*4.0D0**12
984 C Following line for constants currently implemented
985 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
986 aad(i,1)=1.5D0*4.0D0**12
987 c aad(i,1)=0.17D0*5.6D0**12
989 C "Soft" SC-p repulsion
991 C Following line for constants currently implemented
992 c aad(i,1)=0.3D0*4.0D0**6
993 C "Hard" SC-p repulsion
994 bad(i,1)=3.0D0*4.0D0**6
995 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1004 C 8/9/01 Read the SC-p interaction constants from file
1007 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1010 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1011 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1012 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1013 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1017 write (iout,*) "Parameters of SC-p interactions:"
1019 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1020 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1025 C Define the constants of the disulfide bridge
1029 c Old arbitrary potential - commented out.
1034 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1035 c energy surface of diethyl disulfide.
1036 c A. Liwo and U. Kozlowska, 11/24/03
1053 write (iout,'(/a)') "Disulfide bridge parameters:"
1054 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1055 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1056 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1057 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1061 111 write (iout,*) "Error reading bending energy parameters."
1063 112 write (iout,*) "Error reading rotamer energy parameters."
1065 113 write (iout,*) "Error reading torsional energy parameters."
1067 114 write (iout,*) "Error reading double torsional energy parameters."
1070 & "Error reading cumulant (multibody energy) parameters."
1072 116 write (iout,*) "Error reading electrostatic energy parameters."
1074 117 write (iout,*) "Error reading side chain interaction parameters."
1076 118 write (iout,*) "Error reading SCp interaction parameters."
1078 119 write (iout,*) "Error reading SCCOR parameters"
1081 call MPI_Finalize(Ierror)
1088 subroutine getenv_loc(var, val)
1089 character(*) var, val
1092 character(2000) line
1095 open (196,file='env',status='old',readonly,shared)
1097 c write(*,*)'looking for ',var
1098 10 read(196,*,err=11,end=11)line
1099 iread=index(line,var)
1100 c write(*,*)iread,' ',var,' ',line
1101 if (iread.eq.0) go to 10
1102 c write(*,*)'---> ',line
1108 iread=iread+ilen(var)+1
1109 read (line(iread:),*,err=12,end=12) val
1110 c write(*,*)'OK: ',var,' = ',val
1116 #elif (defined CRAY)
1117 integer lennam,lenval,ierror
1119 c getenv using a POSIX call, useful on the T3D
1120 c Sept 1996, comment out error check on advice of H. Pritchard
1123 if(lennam.le.0) stop '--error calling getenv--'
1124 call pxfgetenv(var,lennam,val,lenval,ierror)
1125 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1127 call getenv(var,val)