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)
152 & 'Parameters of the virtual-bond valence angles:'
153 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
154 & ' ATHETA0 ',' A1 ',' A2 ',
157 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
158 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
160 write (iout,'(/a/9x,5a/79(1h-))')
161 & 'Parameters of the expression for sigma(theta_c):',
162 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
163 & ' ALPH3 ',' SIGMA0C '
165 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
166 & (polthet(j,i),j=0,3),sigc0(i)
168 write (iout,'(/a/9x,5a/79(1h-))')
169 & 'Parameters of the second gaussian:',
170 & ' THETA0 ',' SIGMA0 ',' G1 ',
173 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
174 & sig0(i),(gthet(j,i),j=1,3)
178 & 'Parameters of the virtual-bond valence angles:'
179 write (iout,'(/a/9x,5a/79(1h-))')
180 & 'Coefficients of expansion',
181 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
182 & ' b1*10^1 ',' b2*10^1 '
184 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
185 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
186 & (10*bthet(j,i,1,1),j=1,2)
188 write (iout,'(/a/9x,5a/79(1h-))')
189 & 'Parameters of the expression for sigma(theta_c):',
190 & ' alpha0 ',' alph1 ',' alph2 ',
191 & ' alhp3 ',' sigma0c '
193 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
194 & (polthet(j,i),j=0,3),sigc0(i)
196 write (iout,'(/a/9x,5a/79(1h-))')
197 & 'Parameters of the second gaussian:',
198 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
201 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
202 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
208 C Read the parameters of Utheta determined from ab initio surfaces
209 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
211 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
212 & ntheterm3,nsingle,ndouble
213 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
214 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
220 aathet(l,i,j,k)=0.0d0
224 bbthet(m,l,i,j,k)=0.0d0
225 ccthet(m,l,i,j,k)=0.0d0
226 ddthet(m,l,i,j,k)=0.0d0
227 eethet(m,l,i,j,k)=0.0d0
233 ffthet(mm,m,l,i,j,k)=0.0d0
234 ggthet(mm,m,l,i,j,k)=0.0d0
244 read (ithep,'(3a)',end=111,err=111) res1,res2,res3
245 read (ithep,*,end=111,err=111) aa0thet(i,j,k)
246 read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
247 read (ithep,*,end=111,err=111)
248 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
249 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
250 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
251 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
252 read (ithep,*,end=111,err=111)
253 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
254 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
255 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
260 C For dummy ends assign glycine-type coefficients of theta-only terms; the
261 C coefficients of theta-and-gamma-dependent terms are zero.
266 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
267 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
269 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
270 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
273 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
275 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
278 C Control printout of the coefficients of virtual-bond-angle potentials
281 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
285 write (iout,'(//4a)')
286 & 'Type ',onelett(i),onelett(j),onelett(k)
287 write (iout,'(//a,10x,a)') " l","a[l]"
288 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
289 write (iout,'(i2,1pe15.5)')
290 & (l,aathet(l,i,j,k),l=1,ntheterm)
292 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
293 & "b",l,"c",l,"d",l,"e",l
295 write (iout,'(i2,4(1pe15.5))') m,
296 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
297 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
301 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
302 & "f+",l,"f-",l,"g+",l,"g-",l
305 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
306 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
307 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
320 C Read the parameters of the probability distribution/energy expression
321 C of the side chains.
324 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
328 dsc_inv(i)=1.0D0/dsc(i)
339 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
340 & ((blower(k,l,1),l=1,k),k=1,3)
341 censc(1,1,-i)=censc(1,1,i)
342 censc(2,1,-i)=censc(2,1,i)
343 censc(3,1,-i)=-censc(3,1,i)
345 read (irotam,*,end=112,err=112) bsc(j,i)
346 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
347 & ((blower(k,l,j),l=1,k),k=1,3)
348 censc(1,j,-i)=censc(1,j,i)
349 censc(2,j,-i)=censc(2,j,i)
350 censc(3,j,-i)=-censc(3,j,i)
351 C BSC is amplitude of Gaussian
358 akl=akl+blower(k,m,j)*blower(l,m,j)
362 if (((k.eq.3).and.(l.ne.3))
363 & .or.((l.eq.3).and.(k.ne.3))) then
364 gaussc(k,l,j,-i)=-akl
365 gaussc(l,k,j,-i)=-akl
377 write (iout,'(/a)') 'Parameters of side-chain local geometry'
382 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
383 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
384 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
385 & 'log h',(bsc(j,i),j=1,nlobi)
386 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
387 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
389 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
390 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
393 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
394 write (iout,'(a,f10.4,4(16x,f10.4))')
395 & 'Center ',(bsc(j,i),j=1,nlobi)
396 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
405 C Read scrot parameters for potentials determined from all-atom AM1 calculations
406 C added by Urszula Kozlowska 07/11/2007
409 read (irotam,*,end=112,err=112)
411 read (irotam,*,end=112,err=112)
414 read(irotam,*,end=112,err=112) sc_parmin(j,i)
423 C Read torsional parameters in old format
425 read (itorp,*,end=113,err=113) ntortyp,nterm_old
426 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
427 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
432 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
438 write (iout,'(/a/)') 'Torsional constants:'
441 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
442 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
448 C Read torsional parameters
450 read (itorp,*,end=113,err=113) ntortyp
451 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
454 itortyp(i)=-itortyp(-i)
456 c write (iout,*) 'ntortyp',ntortyp
458 do j=-ntortyp+1,ntortyp-1
459 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
461 nterm(-i,-j,iblock)=nterm(i,j,iblock)
462 nlor(-i,-j,iblock)=nlor(i,j,iblock)
465 do k=1,nterm(i,j,iblock)
466 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
468 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
469 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
470 v0ij=v0ij+si*v1(k,i,j,iblock)
472 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
473 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
474 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
476 do k=1,nlor(i,j,iblock)
477 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
478 & vlor2(k,i,j),vlor3(k,i,j)
479 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
482 v0(-i,-j,iblock)=v0ij
488 write (iout,'(/a/)') 'Torsional constants:'
491 write (iout,*) 'ityp',i,' jtyp',j
492 write (iout,*) 'Fourier constants'
493 do k=1,nterm(i,j,iblock)
494 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
497 write (iout,*) 'Lorenz constants'
498 do k=1,nlor(i,j,iblock)
499 write (iout,'(3(1pe15.5))')
500 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
507 C 6/23/01 Read parameters for double torsionals
511 do j=-ntortyp+1,ntortyp-1
512 do k=-ntortyp+1,ntortyp-1
513 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
514 c write (iout,*) "OK onelett",
517 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
518 & .or. t3.ne.toronelet(k)) then
519 write (iout,*) "Error in double torsional parameter file",
522 call MPI_Finalize(Ierror)
524 stop "Error in double torsional parameter file"
526 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
527 & ntermd_2(i,j,k,iblock)
528 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
529 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
530 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
531 & ntermd_1(i,j,k,iblock))
532 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
533 & ntermd_1(i,j,k,iblock))
534 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
535 & ntermd_1(i,j,k,iblock))
536 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
537 & ntermd_1(i,j,k,iblock))
538 C Martix of D parameters for one dimesional foureir series
539 do l=1,ntermd_1(i,j,k,iblock)
540 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
541 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
542 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
543 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
544 c write(iout,*) "whcodze" ,
545 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
547 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
548 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
549 & v2s(m,l,i,j,k,iblock),
550 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
551 C Martix of D parameters for two dimesional fourier series
552 do l=1,ntermd_2(i,j,k,iblock)
554 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
555 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
556 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
557 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
566 write (iout,*) 'Constants for double torsionals'
569 do j=-ntortyp+1,ntortyp-1
570 do k=-ntortyp+1,ntortyp-1
571 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
572 & ' nsingle',ntermd_1(i,j,k,iblock),
573 & ' ndouble',ntermd_2(i,j,k,iblock)
575 write (iout,*) 'Single angles:'
576 do l=1,ntermd_1(i,j,k,iblock)
577 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
578 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
579 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
580 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
583 write (iout,*) 'Pairs of angles:'
584 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
585 do l=1,ntermd_2(i,j,k,iblock)
586 write (iout,'(i5,20f10.5)')
587 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
590 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
591 do l=1,ntermd_2(i,j,k,iblock)
592 write (iout,'(i5,20f10.5)')
593 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
594 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
604 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
605 C correlation energies.
607 read (isccor,*,end=119,err=119) nterm_sccor
612 read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
619 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
622 write (iout,*) 'ityp',i,' jtyp',j
624 write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
630 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
631 C interaction energy of the Gly, Ala, and Pro prototypes.
635 write (iout,*) "Coefficients of the cumulants"
637 read (ifourier,*) nloctyp
639 read (ifourier,*,end=115,err=115)
640 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
642 write (iout,*) 'Type',i
643 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
685 c Ctilde(1,1,i)=0.0d0
686 c Ctilde(1,2,i)=0.0d0
687 c Ctilde(2,1,i)=0.0d0
688 c Ctilde(2,2,i)=0.0d0
710 c Dtilde(1,1,i)=0.0d0
711 c Dtilde(1,2,i)=0.0d0
712 c Dtilde(2,1,i)=0.0d0
713 c Dtilde(2,2,i)=0.0d0
714 EE(1,1,i)= b(10)+b(11)
715 EE(2,2,i)=-b(10)+b(11)
716 EE(2,1,i)= b(12)-b(13)
717 EE(1,2,i)= b(12)+b(13)
718 EE(1,1,-i)= b(10)+b(11)
719 EE(2,2,-i)=-b(10)+b(11)
720 EE(2,1,-i)=-b(12)+b(13)
721 EE(1,2,-i)=-b(12)-b(13)
727 c ee(2,1,i)=ee(1,2,i)
731 write (iout,*) 'Type',i
733 write(iout,*) B1(1,i),B1(2,i)
735 write(iout,*) B2(1,i),B2(2,i)
738 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
742 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
746 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
752 C Read electrostatic-interaction parameters
756 write (iout,'(/a)') 'Electrostatic interaction constants:'
757 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
758 & 'IT','JT','APP','BPP','AEL6','AEL3'
760 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
761 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
762 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
763 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
768 app (i,j)=epp(i,j)*rri*rri
769 bpp (i,j)=-2.0D0*epp(i,j)*rri
770 ael6(i,j)=elpp6(i,j)*4.2D0**6
771 ael3(i,j)=elpp3(i,j)*4.2D0**3
772 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
773 & ael6(i,j),ael3(i,j)
777 C Read side-chain interaction parameters.
779 read (isidep,*,end=117,err=117) ipot,expon
780 if (ipot.lt.1 .or. ipot.gt.5) then
781 write (iout,'(2a)') 'Error while reading SC interaction',
782 & 'potential file - unknown potential type.'
784 call MPI_Finalize(Ierror)
790 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
791 & ', exponents are ',expon,2*expon
792 goto (10,20,30,30,40) ipot
793 C----------------------- LJ potential ---------------------------------
794 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
795 & (sigma0(i),i=1,ntyp)
797 write (iout,'(/a/)') 'Parameters of the LJ potential:'
798 write (iout,'(a/)') 'The epsilon array:'
799 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
800 write (iout,'(/a)') 'One-body parameters:'
801 write (iout,'(a,4x,a)') 'residue','sigma'
802 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
805 C----------------------- LJK potential --------------------------------
806 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
807 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
809 write (iout,'(/a/)') 'Parameters of the LJK potential:'
810 write (iout,'(a/)') 'The epsilon array:'
811 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
812 write (iout,'(/a)') 'One-body parameters:'
813 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
814 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
818 C---------------------- GB or BP potential -----------------------------
819 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
820 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
822 C For the GB potential convert sigma'**2 into chi'
825 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
829 write (iout,'(/a/)') 'Parameters of the BP potential:'
830 write (iout,'(a/)') 'The epsilon array:'
831 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
832 write (iout,'(/a)') 'One-body parameters:'
833 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
835 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
836 & chip(i),alp(i),i=1,ntyp)
839 C--------------------- GBV potential -----------------------------------
840 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
841 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
842 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
844 write (iout,'(/a/)') 'Parameters of the GBV 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,5a)') 'residue',' sigma ',' r0 ',
849 & 's||/s_|_^2',' chip ',' alph '
850 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
851 & sigii(i),chip(i),alp(i),i=1,ntyp)
855 C-----------------------------------------------------------------------
856 C Calculate the "working" parameters of SC interactions.
864 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
865 sigma(j,i)=sigma(i,j)
866 rs0(i,j)=dwa16*sigma(i,j)
870 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
871 & 'Working parameters of the SC interactions:',
872 & ' a ',' b ',' augm ',' sigma ',' r0 ',
877 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
886 sigeps=dsign(1.0D0,epsij)
888 aa(i,j)=epsij*rrij*rrij
889 bb(i,j)=-sigeps*epsij*rrij
897 ratsig1=sigt2sq/sigt1sq
898 ratsig2=1.0D0/ratsig1
899 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
900 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
901 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
905 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
906 sigmaii(i,j)=rsum_max
907 sigmaii(j,i)=rsum_max
909 c sigmaii(i,j)=r0(i,j)
910 c sigmaii(j,i)=r0(i,j)
912 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
913 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
914 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
915 augm(i,j)=epsij*r_augm**(2*expon)
916 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
923 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
924 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
925 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
931 C Define the SC-p interaction constants (hard-coded; old style)
934 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
936 c aad(i,1)=0.3D0*4.0D0**12
937 C Following line for constants currently implemented
938 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
939 aad(i,1)=1.5D0*4.0D0**12
940 c aad(i,1)=0.17D0*5.6D0**12
942 C "Soft" SC-p repulsion
944 C Following line for constants currently implemented
945 c aad(i,1)=0.3D0*4.0D0**6
946 C "Hard" SC-p repulsion
947 bad(i,1)=3.0D0*4.0D0**6
948 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
957 C 8/9/01 Read the SC-p interaction constants from file
960 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
963 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
964 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
965 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
966 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
970 write (iout,*) "Parameters of SC-p interactions:"
972 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
973 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
978 C Define the constants of the disulfide bridge
982 c Old arbitrary potential - commented out.
987 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
988 c energy surface of diethyl disulfide.
989 c A. Liwo and U. Kozlowska, 11/24/03
1006 write (iout,'(/a)') "Disulfide bridge parameters:"
1007 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1008 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1009 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1010 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1014 111 write (iout,*) "Error reading bending energy parameters."
1016 112 write (iout,*) "Error reading rotamer energy parameters."
1018 113 write (iout,*) "Error reading torsional energy parameters."
1020 114 write (iout,*) "Error reading double torsional energy parameters."
1023 & "Error reading cumulant (multibody energy) parameters."
1025 116 write (iout,*) "Error reading electrostatic energy parameters."
1027 117 write (iout,*) "Error reading side chain interaction parameters."
1029 118 write (iout,*) "Error reading SCp interaction parameters."
1031 119 write (iout,*) "Error reading SCCOR parameters"
1034 call MPI_Finalize(Ierror)
1041 subroutine getenv_loc(var, val)
1042 character(*) var, val
1045 character(2000) line
1048 open (196,file='env',status='old',readonly,shared)
1050 c write(*,*)'looking for ',var
1051 10 read(196,*,err=11,end=11)line
1052 iread=index(line,var)
1053 c write(*,*)iread,' ',var,' ',line
1054 if (iread.eq.0) go to 10
1055 c write(*,*)'---> ',line
1061 iread=iread+ilen(var)+1
1062 read (line(iread:),*,err=12,end=12) val
1063 c write(*,*)'OK: ',var,' = ',val
1069 #elif (defined CRAY)
1070 integer lennam,lenval,ierror
1072 c getenv using a POSIX call, useful on the T3D
1073 c Sept 1996, comment out error check on advice of H. Pritchard
1076 if(lennam.le.0) stop '--error calling getenv--'
1077 call pxfgetenv(var,lennam,val,lenval,ierror)
1078 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1080 call getenv(var,val)