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)
288 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)
317 read (irotam,*,end=112,err=112) bsc(j,i)
318 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
319 & ((blower(k,l,j),l=1,k),k=1,3)
326 akl=akl+blower(k,m,j)*blower(l,m,j)
337 write (iout,'(/a)') 'Parameters of side-chain local geometry'
342 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
343 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
344 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
345 & 'log h',(bsc(j,i),j=1,nlobi)
346 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
347 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
349 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
350 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
353 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
354 write (iout,'(a,f10.4,4(16x,f10.4))')
355 & 'Center ',(bsc(j,i),j=1,nlobi)
356 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
365 C Read scrot parameters for potentials determined from all-atom AM1 calculations
366 C added by Urszula Kozlowska 07/11/2007
369 read (irotam,*,end=112,err=112)
371 read (irotam,*,end=112,err=112)
374 read(irotam,*,end=112,err=112) sc_parmin(j,i)
379 C Read the parameters of the probability distribution/energy expression
380 C of the side chains.
383 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
387 dsc_inv(i)=1.0D0/dsc(i)
398 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
399 & ((blower(k,l,1),l=1,k),k=1,3)
401 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
402 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
403 & ((blower(k,l,j),l=1,k),k=1,3)
410 akl=akl+blower(k,m,j)*blower(l,m,j)
425 C Read torsional parameters in old format
427 read (itorp,*,end=113,err=113) ntortyp,nterm_old
428 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
429 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
434 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
440 write (iout,'(/a/)') 'Torsional constants:'
443 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
444 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
450 C Read torsional parameters
452 read (itorp,*,end=113,err=113) ntortyp
453 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
454 c write (iout,*) 'ntortyp',ntortyp
457 read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
461 read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
462 v0ij=v0ij+si*v1(k,i,j)
466 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
467 & vlor2(k,i,j),vlor3(k,i,j)
468 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
475 write (iout,'(/a/)') 'Torsional constants:'
478 write (iout,*) 'ityp',i,' jtyp',j
479 write (iout,*) 'Fourier constants'
481 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
483 write (iout,*) 'Lorenz constants'
485 write (iout,'(3(1pe15.5))')
486 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
492 C 6/23/01 Read parameters for double torsionals
496 do j=-ntortyp+1,ntortyp-1
497 do k=-ntortyp+1,ntortyp-1
498 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
499 write (iout,*) "OK onelett",
502 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
503 & .or. t3.ne.toronelet(k)) then
504 write (iout,*) "Error in double torsional parameter file",
507 call MPI_Finalize(Ierror)
509 stop "Error in double torsional parameter file"
511 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
512 & ntermd_2(i,j,k,iblock)
513 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
514 & ntermd_1(i,j,k,iblock))
515 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
516 & ntermd_1(i,j,k,iblock))
517 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
518 & ntermd_1(i,j,k,iblock))
519 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
520 & ntermd_1(i,j,k,iblock))
521 C Matrix of D parameters for one dimesional foureir series
522 do l=1, ntermd_1(i,j,k,iblock)
523 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
524 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
525 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
526 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
528 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
529 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
530 & v2s(m,l,i,j,k,iblock),
531 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
532 C Matrix of D parameters for two dimesional fourier series
533 do l=1,ntermd_2(i,j,k,iblock)
535 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
536 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
537 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
538 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
547 write (iout,*) 'Constants for double torsionals'
549 do j=-ntortyp,ntortyp
550 do k=-ntortyp,ntortyp
551 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
552 & ' nsingle',ntermd_1(i,j,k,iblock),
553 & ' ndouble',ntermd_2(i,j,k,iblock)
555 write (iout,*) 'Single angles:'
556 do l=1,ntermd_1(i,j,k,iblock)
557 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
558 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
559 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
562 write (iout,*) 'Pairs of angles:'
563 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
564 do l=1,ntermd_2(i,j,k,iblock)
565 write (iout,'(i5,20f10.5)')
566 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
569 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
570 do l=1,ntermd_2(i,j,k,iblock)
571 write (iout,'(i5,20f10.5)')
572 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
573 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
582 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
583 C correlation energies.
585 read (isccor,*,end=119,err=119) nterm_sccor
590 read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
597 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
600 write (iout,*) 'ityp',i,' jtyp',j
602 write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
608 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
609 C interaction energy of the Gly, Ala, and Pro prototypes.
613 write (iout,*) "Coefficients of the cumulants"
615 read (ifourier,*) nloctyp
618 read (ifourier,*,end=115,err=115)
619 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
621 write (iout,*) 'Type',i
622 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
648 c Ctilde(1,1,i)=0.0d0
649 c Ctilde(1,2,i)=0.0d0
650 c Ctilde(2,1,i)=0.0d0
651 c Ctilde(2,2,i)=0.0d0
664 c Dtilde(1,1,i)=0.0d0
665 c Dtilde(1,2,i)=0.0d0
666 c Dtilde(2,1,i)=0.0d0
667 c Dtilde(2,2,i)=0.0d0
668 EE(1,1,i)= b(10)+b(11)
669 EE(2,2,i)=-b(10)+b(11)
670 EE(2,1,i)= b(12)-b(13)
671 EE(1,2,i)= b(12)+b(13)
676 c ee(2,1,i)=ee(1,2,i)
680 write (iout,*) 'Type',i
682 write(iout,*) B1(1,i),B1(2,i)
684 write(iout,*) B2(1,i),B2(2,i)
687 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
691 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
695 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
700 C Read electrostatic-interaction parameters
704 write (iout,'(/a)') 'Electrostatic interaction constants:'
705 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
706 & 'IT','JT','APP','BPP','AEL6','AEL3'
708 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
709 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
710 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
711 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
716 app (i,j)=epp(i,j)*rri*rri
717 bpp (i,j)=-2.0D0*epp(i,j)*rri
718 ael6(i,j)=elpp6(i,j)*4.2D0**6
719 ael3(i,j)=elpp3(i,j)*4.2D0**3
720 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
721 & ael6(i,j),ael3(i,j)
725 C Read side-chain interaction parameters.
727 read (isidep,*,end=117,err=117) ipot,expon
728 if (ipot.lt.1 .or. ipot.gt.5) then
729 write (iout,'(2a)') 'Error while reading SC interaction',
730 & 'potential file - unknown potential type.'
732 call MPI_Finalize(Ierror)
738 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
739 & ', exponents are ',expon,2*expon
740 goto (10,20,30,30,40) ipot
741 C----------------------- LJ potential ---------------------------------
742 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
743 & (sigma0(i),i=1,ntyp)
745 write (iout,'(/a/)') 'Parameters of the LJ potential:'
746 write (iout,'(a/)') 'The epsilon array:'
747 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
748 write (iout,'(/a)') 'One-body parameters:'
749 write (iout,'(a,4x,a)') 'residue','sigma'
750 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
753 C----------------------- LJK potential --------------------------------
754 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
755 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
757 write (iout,'(/a/)') 'Parameters of the LJK potential:'
758 write (iout,'(a/)') 'The epsilon array:'
759 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
760 write (iout,'(/a)') 'One-body parameters:'
761 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
762 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
766 C---------------------- GB or BP potential -----------------------------
767 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
768 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
770 C For the GB potential convert sigma'**2 into chi'
773 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
777 write (iout,'(/a/)') 'Parameters of the BP potential:'
778 write (iout,'(a/)') 'The epsilon array:'
779 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
780 write (iout,'(/a)') 'One-body parameters:'
781 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
783 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
784 & chip(i),alp(i),i=1,ntyp)
787 C--------------------- GBV potential -----------------------------------
788 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
789 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
790 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
792 write (iout,'(/a/)') 'Parameters of the GBV potential:'
793 write (iout,'(a/)') 'The epsilon array:'
794 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
795 write (iout,'(/a)') 'One-body parameters:'
796 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
797 & 's||/s_|_^2',' chip ',' alph '
798 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
799 & sigii(i),chip(i),alp(i),i=1,ntyp)
803 C-----------------------------------------------------------------------
804 C Calculate the "working" parameters of SC interactions.
812 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
813 sigma(j,i)=sigma(i,j)
814 rs0(i,j)=dwa16*sigma(i,j)
818 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
819 & 'Working parameters of the SC interactions:',
820 & ' a ',' b ',' augm ',' sigma ',' r0 ',
825 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
834 sigeps=dsign(1.0D0,epsij)
836 aa(i,j)=epsij*rrij*rrij
837 bb(i,j)=-sigeps*epsij*rrij
845 ratsig1=sigt2sq/sigt1sq
846 ratsig2=1.0D0/ratsig1
847 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
848 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
849 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
853 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
854 sigmaii(i,j)=rsum_max
855 sigmaii(j,i)=rsum_max
857 c sigmaii(i,j)=r0(i,j)
858 c sigmaii(j,i)=r0(i,j)
860 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
861 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
862 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
863 augm(i,j)=epsij*r_augm**(2*expon)
864 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
871 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
872 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
873 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
879 C Define the SC-p interaction constants (hard-coded; old style)
882 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
884 c aad(i,1)=0.3D0*4.0D0**12
885 C Following line for constants currently implemented
886 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
887 aad(i,1)=1.5D0*4.0D0**12
888 c aad(i,1)=0.17D0*5.6D0**12
890 C "Soft" SC-p repulsion
892 C Following line for constants currently implemented
893 c aad(i,1)=0.3D0*4.0D0**6
894 C "Hard" SC-p repulsion
895 bad(i,1)=3.0D0*4.0D0**6
896 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
905 C 8/9/01 Read the SC-p interaction constants from file
908 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
911 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
912 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
913 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
914 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
918 write (iout,*) "Parameters of SC-p interactions:"
920 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
921 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
926 C Define the constants of the disulfide bridge
930 c Old arbitrary potential - commented out.
935 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
936 c energy surface of diethyl disulfide.
937 c A. Liwo and U. Kozlowska, 11/24/03
954 write (iout,'(/a)') "Disulfide bridge parameters:"
955 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
956 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
957 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
958 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
962 111 write (iout,*) "Error reading bending energy parameters."
964 112 write (iout,*) "Error reading rotamer energy parameters."
966 113 write (iout,*) "Error reading torsional energy parameters."
968 114 write (iout,*) "Error reading double torsional energy parameters."
971 & "Error reading cumulant (multibody energy) parameters."
973 116 write (iout,*) "Error reading electrostatic energy parameters."
975 117 write (iout,*) "Error reading side chain interaction parameters."
977 118 write (iout,*) "Error reading SCp interaction parameters."
979 119 write (iout,*) "Error reading SCCOR parameters"
982 call MPI_Finalize(Ierror)
989 subroutine getenv_loc(var, val)
990 character(*) var, val
996 open (196,file='env',status='old',readonly,shared)
998 c write(*,*)'looking for ',var
999 10 read(196,*,err=11,end=11)line
1000 iread=index(line,var)
1001 c write(*,*)iread,' ',var,' ',line
1002 if (iread.eq.0) go to 10
1003 c write(*,*)'---> ',line
1009 iread=iread+ilen(var)+1
1010 read (line(iread:),*,err=12,end=12) val
1011 c write(*,*)'OK: ',var,' = ',val
1017 #elif (defined CRAY)
1018 integer lennam,lenval,ierror
1020 c getenv using a POSIX call, useful on the T3D
1021 c Sept 1996, comment out error check on advice of H. Pritchard
1024 if(lennam.le.0) stop '--error calling getenv--'
1025 call pxfgetenv(var,lennam,val,lenval,ierror)
1026 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1028 call getenv(var,val)