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 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))
603 C Read of Side-chain backbone correlation parameters
604 C Modified 11 May 2012 by Adasko
607 read (isccor,*,end=119,err=119) nsccortyp
609 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
611 isccortyp(i)=-isccortyp(-i)
613 iscprol=isccortyp(20)
614 c write (iout,*) 'ntortyp',ntortyp
616 cc maxinter is maximum interaction sites
620 read (isccor,*,end=119,err=119)
621 &nterm_sccor(i,j),nlor_sccor(i,j)
627 nterm_sccor(-i,j)=nterm_sccor(i,j)
628 nterm_sccor(-i,-j)=nterm_sccor(i,j)
629 nterm_sccor(i,-j)=nterm_sccor(i,j)
630 do k=1,nterm_sccor(i,j)
631 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
633 if (j.eq.iscprol) then
634 if (i.eq.isccortyp(10)) then
635 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
636 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
638 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
639 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
640 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
641 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
642 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
643 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
644 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
645 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
648 if (i.eq.isccortyp(10)) then
649 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
650 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
652 if (j.eq.isccortyp(10)) then
653 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
654 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
656 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
657 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
658 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
659 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
660 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
661 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
665 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
666 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
667 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
668 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
671 do k=1,nlor_sccor(i,j)
672 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
673 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
674 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
675 &(1+vlor3sccor(k,i,j)**2)
677 v0sccor(l,i,j)=v0ijsccor
678 v0sccor(l,-i,j)=v0ijsccor1
679 v0sccor(l,i,-j)=v0ijsccor2
680 v0sccor(l,-i,-j)=v0ijsccor3
686 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
687 c write (iout,*) 'ntortyp',ntortyp
689 cc maxinter is maximum interaction sites
693 read (isccor,*,end=113,err=113)
694 & nterm_sccor(i,j),nlor_sccor(i,j)
698 do k=1,nterm_sccor(i,j)
699 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
701 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
704 do k=1,nlor_sccor(i,j)
705 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
706 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
707 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
708 &(1+vlor3sccor(k,i,j)**2)
710 v0sccor(i,j)=v0ijsccor
718 write (iout,'(/a/)') 'Torsional constants:'
721 write (iout,*) 'ityp',i,' jtyp',j
722 write (iout,*) 'Fourier constants'
723 do k=1,nterm_sccor(i,j)
724 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
726 write (iout,*) 'Lorenz constants'
727 do k=1,nlor_sccor(i,j)
728 write (iout,'(3(1pe15.5))')
729 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
736 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
737 C interaction energy of the Gly, Ala, and Pro prototypes.
741 write (iout,*) "Coefficients of the cumulants"
743 read (ifourier,*) nloctyp
745 read (ifourier,*,end=115,err=115)
746 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
748 write (iout,*) 'Type',i
749 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
791 c Ctilde(1,1,i)=0.0d0
792 c Ctilde(1,2,i)=0.0d0
793 c Ctilde(2,1,i)=0.0d0
794 c Ctilde(2,2,i)=0.0d0
816 c Dtilde(1,1,i)=0.0d0
817 c Dtilde(1,2,i)=0.0d0
818 c Dtilde(2,1,i)=0.0d0
819 c Dtilde(2,2,i)=0.0d0
820 EE(1,1,i)= b(10)+b(11)
821 EE(2,2,i)=-b(10)+b(11)
822 EE(2,1,i)= b(12)-b(13)
823 EE(1,2,i)= b(12)+b(13)
824 EE(1,1,-i)= b(10)+b(11)
825 EE(2,2,-i)=-b(10)+b(11)
826 EE(2,1,-i)=-b(12)+b(13)
827 EE(1,2,-i)=-b(12)-b(13)
833 c ee(2,1,i)=ee(1,2,i)
837 write (iout,*) 'Type',i
839 write(iout,*) B1(1,i),B1(2,i)
841 write(iout,*) B2(1,i),B2(2,i)
844 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
848 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
852 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
858 C Read electrostatic-interaction parameters
862 write (iout,'(/a)') 'Electrostatic interaction constants:'
863 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
864 & 'IT','JT','APP','BPP','AEL6','AEL3'
866 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
867 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
868 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
869 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
874 app (i,j)=epp(i,j)*rri*rri
875 bpp (i,j)=-2.0D0*epp(i,j)*rri
876 ael6(i,j)=elpp6(i,j)*4.2D0**6
877 ael3(i,j)=elpp3(i,j)*4.2D0**3
878 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
879 & ael6(i,j),ael3(i,j)
883 C Read side-chain interaction parameters.
885 read (isidep,*,end=117,err=117) ipot,expon
886 if (ipot.lt.1 .or. ipot.gt.5) then
887 write (iout,'(2a)') 'Error while reading SC interaction',
888 & 'potential file - unknown potential type.'
890 call MPI_Finalize(Ierror)
896 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
897 & ', exponents are ',expon,2*expon
898 goto (10,20,30,30,40) ipot
899 C----------------------- LJ potential ---------------------------------
900 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
901 & (sigma0(i),i=1,ntyp)
903 write (iout,'(/a/)') 'Parameters of the LJ potential:'
904 write (iout,'(a/)') 'The epsilon array:'
905 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
906 write (iout,'(/a)') 'One-body parameters:'
907 write (iout,'(a,4x,a)') 'residue','sigma'
908 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
911 C----------------------- LJK potential --------------------------------
912 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
913 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
915 write (iout,'(/a/)') 'Parameters of the LJK potential:'
916 write (iout,'(a/)') 'The epsilon array:'
917 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
918 write (iout,'(/a)') 'One-body parameters:'
919 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
920 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
924 C---------------------- GB or BP potential -----------------------------
926 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
928 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
929 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
930 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
931 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
932 C For the GB potential convert sigma'**2 into chi'
935 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
939 write (iout,'(/a/)') 'Parameters of the BP 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,4a)') 'residue',' sigma ','s||/s_|_^2',
945 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
946 & chip(i),alp(i),i=1,ntyp)
949 C--------------------- GBV potential -----------------------------------
950 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
951 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
952 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
954 write (iout,'(/a/)') 'Parameters of the GBV potential:'
955 write (iout,'(a/)') 'The epsilon array:'
956 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
957 write (iout,'(/a)') 'One-body parameters:'
958 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
959 & 's||/s_|_^2',' chip ',' alph '
960 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
961 & sigii(i),chip(i),alp(i),i=1,ntyp)
965 C-----------------------------------------------------------------------
966 C Calculate the "working" parameters of SC interactions.
974 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
975 sigma(j,i)=sigma(i,j)
976 rs0(i,j)=dwa16*sigma(i,j)
980 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
981 & 'Working parameters of the SC interactions:',
982 & ' a ',' b ',' augm ',' sigma ',' r0 ',
987 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
996 sigeps=dsign(1.0D0,epsij)
998 aa(i,j)=epsij*rrij*rrij
999 bb(i,j)=-sigeps*epsij*rrij
1003 sigt1sq=sigma0(i)**2
1004 sigt2sq=sigma0(j)**2
1007 ratsig1=sigt2sq/sigt1sq
1008 ratsig2=1.0D0/ratsig1
1009 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1010 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1011 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1015 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1016 sigmaii(i,j)=rsum_max
1017 sigmaii(j,i)=rsum_max
1019 c sigmaii(i,j)=r0(i,j)
1020 c sigmaii(j,i)=r0(i,j)
1022 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1023 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1024 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1025 augm(i,j)=epsij*r_augm**(2*expon)
1026 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1033 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1034 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1035 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1041 C Define the SC-p interaction constants (hard-coded; old style)
1044 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1046 c aad(i,1)=0.3D0*4.0D0**12
1047 C Following line for constants currently implemented
1048 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1049 aad(i,1)=1.5D0*4.0D0**12
1050 c aad(i,1)=0.17D0*5.6D0**12
1052 C "Soft" SC-p repulsion
1054 C Following line for constants currently implemented
1055 c aad(i,1)=0.3D0*4.0D0**6
1056 C "Hard" SC-p repulsion
1057 bad(i,1)=3.0D0*4.0D0**6
1058 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1067 C 8/9/01 Read the SC-p interaction constants from file
1070 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1073 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1074 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1075 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1076 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1080 write (iout,*) "Parameters of SC-p interactions:"
1082 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1083 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1088 C Define the constants of the disulfide bridge
1092 c Old arbitrary potential - commented out.
1097 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1098 c energy surface of diethyl disulfide.
1099 c A. Liwo and U. Kozlowska, 11/24/03
1116 write (iout,'(/a)') "Disulfide bridge parameters:"
1117 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1118 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1119 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1120 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1124 111 write (iout,*) "Error reading bending energy parameters."
1126 112 write (iout,*) "Error reading rotamer energy parameters."
1128 113 write (iout,*) "Error reading torsional energy parameters."
1130 114 write (iout,*) "Error reading double torsional energy parameters."
1133 & "Error reading cumulant (multibody energy) parameters."
1135 116 write (iout,*) "Error reading electrostatic energy parameters."
1137 117 write (iout,*) "Error reading side chain interaction parameters."
1139 118 write (iout,*) "Error reading SCp interaction parameters."
1141 119 write (iout,*) "Error reading SCCOR parameters"
1144 call MPI_Finalize(Ierror)
1151 subroutine getenv_loc(var, val)
1152 character(*) var, val
1155 character(2000) line
1158 open (196,file='env',status='old',readonly,shared)
1160 c write(*,*)'looking for ',var
1161 10 read(196,*,err=11,end=11)line
1162 iread=index(line,var)
1163 c write(*,*)iread,' ',var,' ',line
1164 if (iread.eq.0) go to 10
1165 c write(*,*)'---> ',line
1171 iread=iread+ilen(var)+1
1172 read (line(iread:),*,err=12,end=12) val
1173 c write(*,*)'OK: ',var,' = ',val
1179 #elif (defined CRAY)
1180 integer lennam,lenval,ierror
1182 c getenv using a POSIX call, useful on the T3D
1183 c Sept 1996, comment out error check on advice of H. Pritchard
1186 if(lennam.le.0) stop '--error calling getenv--'
1187 call pxfgetenv(var,lennam,val,lenval,ierror)
1188 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1190 call getenv(var,val)