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),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)
117 & 'Parameters of the virtual-bond valence angles:'
118 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
119 & ' ATHETA0 ',' A1 ',' A2 ',
122 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
123 & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
125 write (iout,'(/a/9x,5a/79(1h-))')
126 & 'Parameters of the expression for sigma(theta_c):',
127 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
128 & ' ALPH3 ',' SIGMA0C '
130 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
131 & (polthet(j,i),j=0,3),sigc0(i)
133 write (iout,'(/a/9x,5a/79(1h-))')
134 & 'Parameters of the second gaussian:',
135 & ' THETA0 ',' SIGMA0 ',' G1 ',
138 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
139 & sig0(i),(gthet(j,i),j=1,3)
143 & 'Parameters of the virtual-bond valence angles:'
144 write (iout,'(/a/9x,5a/79(1h-))')
145 & 'Coefficients of expansion',
146 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
147 & ' b1*10^1 ',' b2*10^1 '
149 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
150 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
152 write (iout,'(/a/9x,5a/79(1h-))')
153 & 'Parameters of the expression for sigma(theta_c):',
154 & ' alpha0 ',' alph1 ',' alph2 ',
155 & ' alhp3 ',' sigma0c '
157 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
158 & (polthet(j,i),j=0,3),sigc0(i)
160 write (iout,'(/a/9x,5a/79(1h-))')
161 & 'Parameters of the second gaussian:',
162 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
165 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
166 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
172 C Read the parameters of Utheta determined from ab initio surfaces
173 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
175 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
176 & ntheterm3,nsingle,ndouble
177 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
178 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
184 aathet(l,i,j,k)=0.0d0
188 bbthet(m,l,i,j,k)=0.0d0
189 ccthet(m,l,i,j,k)=0.0d0
190 ddthet(m,l,i,j,k)=0.0d0
191 eethet(m,l,i,j,k)=0.0d0
197 ffthet(mm,m,l,i,j,k)=0.0d0
198 ggthet(mm,m,l,i,j,k)=0.0d0
208 read (ithep,'(3a)',end=111,err=111) res1,res2,res3
209 read (ithep,*,end=111,err=111) aa0thet(i,j,k)
210 read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
211 read (ithep,*,end=111,err=111)
212 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
213 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
214 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
215 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
216 read (ithep,*,end=111,err=111)
217 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
218 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
219 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
224 C For dummy ends assign glycine-type coefficients of theta-only terms; the
225 C coefficients of theta-and-gamma-dependent terms are zero.
230 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
231 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
233 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
234 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
237 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
239 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
242 C Control printout of the coefficients of virtual-bond-angle potentials
245 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
249 write (iout,'(//4a)')
250 & 'Type ',onelett(i),onelett(j),onelett(k)
251 write (iout,'(//a,10x,a)') " l","a[l]"
252 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
253 write (iout,'(i2,1pe15.5)')
254 & (l,aathet(l,i,j,k),l=1,ntheterm)
256 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
257 & "b",l,"c",l,"d",l,"e",l
259 write (iout,'(i2,4(1pe15.5))') m,
260 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
261 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
265 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
266 & "f+",l,"f-",l,"g+",l,"g-",l
269 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
270 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
271 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
280 write (2,*) "Start reading THETA_PDB"
282 read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
284 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
285 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
286 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
289 write (2,*) "End reading THETA_PDB"
295 C Read the parameters of the probability distribution/energy expression
296 C of the side chains.
299 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
303 dsc_inv(i)=1.0D0/dsc(i)
314 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
315 & ((blower(k,l,1),l=1,k),k=1,3)
316 censc(1,1,-i)=censc(1,1,i)
317 censc(2,1,-i)=censc(2,1,i)
318 censc(3,1,-i)=-censc(3,1,i)
321 read (irotam,*,end=112,err=112) bsc(j,i)
322 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
323 & ((blower(k,l,j),l=1,k),k=1,3)
324 censc(1,j,-i)=censc(1,j,i)
325 censc(2,j,-i)=censc(2,j,i)
326 censc(3,j,-i)=-censc(3,j,i)
327 C BSC is amplitude of Gaussian
334 akl=akl+blower(k,m,j)*blower(l,m,j)
338 if (((k.eq.3).and.(l.ne.3))
339 & .or.((l.eq.3).and.(k.ne.3))) then
340 gaussc(k,l,j,-i)=-akl
341 gaussc(l,k,j,-i)=-akl
353 write (iout,'(/a)') 'Parameters of side-chain local geometry'
358 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
359 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
360 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
361 & 'log h',(bsc(j,i),j=1,nlobi)
362 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
363 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
365 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
366 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
369 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
370 write (iout,'(a,f10.4,4(16x,f10.4))')
371 & 'Center ',(bsc(j,i),j=1,nlobi)
372 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
381 C Read scrot parameters for potentials determined from all-atom AM1 calculations
382 C added by Urszula Kozlowska 07/11/2007
385 read (irotam,*,end=112,err=112)
387 read (irotam,*,end=112,err=112)
390 read(irotam,*,end=112,err=112) sc_parmin(j,i)
395 C Read the parameters of the probability distribution/energy expression
396 C of the side chains.
399 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
403 dsc_inv(i)=1.0D0/dsc(i)
414 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
415 & ((blower(k,l,1),l=1,k),k=1,3)
417 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
418 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
419 & ((blower(k,l,j),l=1,k),k=1,3)
426 akl=akl+blower(k,m,j)*blower(l,m,j)
441 C Read torsional parameters in old format
443 read (itorp,*,end=113,err=113) ntortyp,nterm_old
444 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
445 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
450 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
456 write (iout,'(/a/)') 'Torsional constants:'
459 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
460 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
466 C Read torsional parameters
468 read (itorp,*,end=113,err=113) ntortyp
469 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
472 itortyp(i)=-itortyp(-i)
474 c write (iout,*) 'ntortyp',ntortyp
476 do j=-ntortyp+1,ntortyp-1
477 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
479 nterm(-i,-j,iblock)=nterm(i,j,iblock)
480 nlor(-i,-j,iblock)=nlor(i,j,iblock)
483 do k=1,nterm(i,j,iblock)
484 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
486 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
487 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
488 v0ij=v0ij+si*v1(k,i,j,iblock)
490 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock),v1(k,-i,-j,iblock)
492 do k=1,nlor(i,j,iblock)
493 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
494 & vlor2(k,i,j),vlor3(k,i,j)
495 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
498 v0(-i,-j,iblock)=v0ij
504 write (iout,'(/a/)') 'Torsional constants:'
507 write (iout,*) 'ityp',i,' jtyp',j
508 write (iout,*) 'Fourier constants'
509 do k=1,nterm(i,j,iblock)
510 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
513 write (iout,*) 'Lorenz constants'
514 do k=1,nlor(i,j,iblock)
515 write (iout,'(3(1pe15.5))')
516 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
522 C 6/23/01 Read parameters for double torsionals
526 do j=-ntortyp+1,ntortyp-1
527 do k=-ntortyp+1,ntortyp-1
528 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
529 write (iout,*) "OK onelett",
532 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
533 & .or. t3.ne.toronelet(k)) then
534 write (iout,*) "Error in double torsional parameter file",
537 call MPI_Finalize(Ierror)
539 stop "Error in double torsional parameter file"
541 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
542 & ntermd_2(i,j,k,iblock)
543 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
544 & ntermd_1(i,j,k,iblock))
545 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
546 & ntermd_1(i,j,k,iblock))
547 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
548 & ntermd_1(i,j,k,iblock))
549 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
550 & ntermd_1(i,j,k,iblock))
551 C Martix of D parameters for one dimesional foureir series
552 do l=1, ntermd_1(i,j,k,iblock)
553 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
554 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
555 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
556 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
558 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
559 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
560 & v2s(m,l,i,j,k,iblock),
561 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
562 C Martix of D parameters for two dimesional fourier series
563 do l=1,ntermd_2(i,j,k,iblock)
565 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
566 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
567 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
568 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
577 write (iout,*) 'Constants for double torsionals'
579 do j=-ntortyp,ntortyp
580 do k=-ntortyp,ntortyp
581 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
582 & ' nsingle',ntermd_1(i,j,k,iblock),
583 & ' ndouble',ntermd_2(i,j,k,iblock)
585 write (iout,*) 'Single angles:'
586 do l=1,ntermd_1(i,j,k,iblock)
587 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
588 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
589 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
590 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
593 write (iout,*) 'Pairs of angles:'
594 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
595 do l=1,ntermd_2(i,j,k,iblock)
596 write (iout,'(i5,20f10.5)')
597 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
600 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
601 do l=1,ntermd_2(i,j,k,iblock)
602 write (iout,'(i5,20f10.5)')
603 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
604 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
612 C Read of Side-chain backbone correlation parameters
613 C Modified 11 May 2012 by Adasko
616 read (isccor,*,end=113,err=113) nsccortyp
618 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
620 read (isccor,*,end=113,err=113) (isccortyp(i),i=-ntyp,ntyp)
622 c write (iout,*) 'ntortyp',ntortyp
624 cc maxinter is maximum interaction sites
628 read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
632 do k=1,nterm_sccor(i,j)
633 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
635 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
638 do k=1,nlor_sccor(i,j)
639 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
640 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
641 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
642 &(1+vlor3sccor(k,i,j)**2)
644 v0sccor(i,j)=v0ijsccor
651 write (iout,'(/a/)') 'Torsional constants:'
654 write (iout,*) 'ityp',i,' jtyp',j
655 write (iout,*) 'Fourier constants'
656 do k=1,nterm_sccor(i,j)
657 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
659 write (iout,*) 'Lorenz constants'
660 do k=1,nlor_sccor(i,j)
661 write (iout,'(3(1pe15.5))')
662 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
669 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
670 C interaction energy of the Gly, Ala, and Pro prototypes.
674 write (iout,*) "Coefficients of the cumulants"
676 read (ifourier,*) nloctyp
678 read (ifourier,*,end=115,err=115)
679 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
681 write (iout,*) 'Type',i
682 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
724 c Ctilde(1,1,i)=0.0d0
725 c Ctilde(1,2,i)=0.0d0
726 c Ctilde(2,1,i)=0.0d0
727 c Ctilde(2,2,i)=0.0d0
749 c Dtilde(1,1,i)=0.0d0
750 c Dtilde(1,2,i)=0.0d0
751 c Dtilde(2,1,i)=0.0d0
752 c Dtilde(2,2,i)=0.0d0
753 EE(1,1,i)= b(10)+b(11)
754 EE(2,2,i)=-b(10)+b(11)
755 EE(2,1,i)= b(12)-b(13)
756 EE(1,2,i)= b(12)+b(13)
757 EE(1,1,-i)= b(10)+b(11)
758 EE(2,2,-i)=-b(10)+b(11)
759 EE(2,1,-i)=-b(12)+b(13)
760 EE(1,2,-i)=-b(12)-b(13)
766 c ee(2,1,i)=ee(1,2,i)
770 write (iout,*) 'Type',i
772 write(iout,*) B1(1,i),B1(2,i)
774 write(iout,*) B2(1,i),B2(2,i)
777 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
781 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
785 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
790 C Read electrostatic-interaction parameters
794 write (iout,'(/a)') 'Electrostatic interaction constants:'
795 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
796 & 'IT','JT','APP','BPP','AEL6','AEL3'
798 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
799 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
800 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
801 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
806 app (i,j)=epp(i,j)*rri*rri
807 bpp (i,j)=-2.0D0*epp(i,j)*rri
808 ael6(i,j)=elpp6(i,j)*4.2D0**6
809 ael3(i,j)=elpp3(i,j)*4.2D0**3
810 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
811 & ael6(i,j),ael3(i,j)
815 C Read side-chain interaction parameters.
817 read (isidep,*,end=117,err=117) ipot,expon
818 if (ipot.lt.1 .or. ipot.gt.5) then
819 write (iout,'(2a)') 'Error while reading SC interaction',
820 & 'potential file - unknown potential type.'
822 call MPI_Finalize(Ierror)
828 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
829 & ', exponents are ',expon,2*expon
830 goto (10,20,30,30,40) ipot
831 C----------------------- LJ potential ---------------------------------
832 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
833 & (sigma0(i),i=1,ntyp)
835 write (iout,'(/a/)') 'Parameters of the LJ potential:'
836 write (iout,'(a/)') 'The epsilon array:'
837 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
838 write (iout,'(/a)') 'One-body parameters:'
839 write (iout,'(a,4x,a)') 'residue','sigma'
840 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
843 C----------------------- LJK potential --------------------------------
844 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
845 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
847 write (iout,'(/a/)') 'Parameters of the LJK potential:'
848 write (iout,'(a/)') 'The epsilon array:'
849 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
850 write (iout,'(/a)') 'One-body parameters:'
851 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
852 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
856 C---------------------- GB or BP potential -----------------------------
857 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
858 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
860 C For the GB potential convert sigma'**2 into chi'
863 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
867 write (iout,'(/a/)') 'Parameters of the BP potential:'
868 write (iout,'(a/)') 'The epsilon array:'
869 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
870 write (iout,'(/a)') 'One-body parameters:'
871 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
873 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
874 & chip(i),alp(i),i=1,ntyp)
877 C--------------------- GBV potential -----------------------------------
878 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
879 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
880 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
882 write (iout,'(/a/)') 'Parameters of the GBV potential:'
883 write (iout,'(a/)') 'The epsilon array:'
884 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
885 write (iout,'(/a)') 'One-body parameters:'
886 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
887 & 's||/s_|_^2',' chip ',' alph '
888 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
889 & sigii(i),chip(i),alp(i),i=1,ntyp)
893 C-----------------------------------------------------------------------
894 C Calculate the "working" parameters of SC interactions.
902 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
903 sigma(j,i)=sigma(i,j)
904 rs0(i,j)=dwa16*sigma(i,j)
908 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
909 & 'Working parameters of the SC interactions:',
910 & ' a ',' b ',' augm ',' sigma ',' r0 ',
915 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
924 sigeps=dsign(1.0D0,epsij)
926 aa(i,j)=epsij*rrij*rrij
927 bb(i,j)=-sigeps*epsij*rrij
935 ratsig1=sigt2sq/sigt1sq
936 ratsig2=1.0D0/ratsig1
937 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
938 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
939 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
943 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
944 sigmaii(i,j)=rsum_max
945 sigmaii(j,i)=rsum_max
947 c sigmaii(i,j)=r0(i,j)
948 c sigmaii(j,i)=r0(i,j)
950 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
951 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
952 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
953 augm(i,j)=epsij*r_augm**(2*expon)
954 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
961 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
962 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
963 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
969 C Define the SC-p interaction constants (hard-coded; old style)
972 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
974 c aad(i,1)=0.3D0*4.0D0**12
975 C Following line for constants currently implemented
976 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
977 aad(i,1)=1.5D0*4.0D0**12
978 c aad(i,1)=0.17D0*5.6D0**12
980 C "Soft" SC-p repulsion
982 C Following line for constants currently implemented
983 c aad(i,1)=0.3D0*4.0D0**6
984 C "Hard" SC-p repulsion
985 bad(i,1)=3.0D0*4.0D0**6
986 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
995 C 8/9/01 Read the SC-p interaction constants from file
998 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1001 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1002 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1003 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1004 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1008 write (iout,*) "Parameters of SC-p interactions:"
1010 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1011 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1016 C Define the constants of the disulfide bridge
1020 c Old arbitrary potential - commented out.
1025 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1026 c energy surface of diethyl disulfide.
1027 c A. Liwo and U. Kozlowska, 11/24/03
1044 write (iout,'(/a)') "Disulfide bridge parameters:"
1045 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1046 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1047 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1048 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1052 111 write (iout,*) "Error reading bending energy parameters."
1054 112 write (iout,*) "Error reading rotamer energy parameters."
1056 113 write (iout,*) "Error reading torsional energy parameters."
1058 114 write (iout,*) "Error reading double torsional energy parameters."
1061 & "Error reading cumulant (multibody energy) parameters."
1063 116 write (iout,*) "Error reading electrostatic energy parameters."
1065 117 write (iout,*) "Error reading side chain interaction parameters."
1067 118 write (iout,*) "Error reading SCp interaction parameters."
1069 119 write (iout,*) "Error reading SCCOR parameters"
1072 call MPI_Finalize(Ierror)
1079 subroutine getenv_loc(var, val)
1080 character(*) var, val
1083 character(2000) line
1086 open (196,file='env',status='old',readonly,shared)
1088 c write(*,*)'looking for ',var
1089 10 read(196,*,err=11,end=11)line
1090 iread=index(line,var)
1091 c write(*,*)iread,' ',var,' ',line
1092 if (iread.eq.0) go to 10
1093 c write(*,*)'---> ',line
1099 iread=iread+ilen(var)+1
1100 read (line(iread:),*,err=12,end=12) val
1101 c write(*,*)'OK: ',var,' = ',val
1107 #elif (defined CRAY)
1108 integer lennam,lenval,ierror
1110 c getenv using a POSIX call, useful on the T3D
1111 c Sept 1996, comment out error check on advice of H. Pritchard
1114 if(lennam.le.0) stop '--error calling getenv--'
1115 call pxfgetenv(var,lennam,val,lenval,ierror)
1116 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1118 call getenv(var,val)