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))
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)
624 nterm_sccor(-i,j)=nterm_sccor(i,j)
625 nterm_sccor(-i,-j)=nterm_sccor(i,j)
626 nterm_sccor(i,-j)=nterm_sccor(i,j)
627 do k=1,nterm_sccor(i,j)
628 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
630 if (j.eq.iscprol) then
631 if (i.eq.isccortyp(10)) then
632 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
633 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
635 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
636 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
637 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
638 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
639 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
640 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
641 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
642 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
645 if (i.eq.isccortyp(10)) then
646 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
647 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
649 if (j.eq.isccortyp(10)) then
650 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
651 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
653 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
654 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
655 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
656 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
657 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
658 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
662 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
665 do k=1,nlor_sccor(i,j)
666 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
667 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
668 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
669 &(1+vlor3sccor(k,i,j)**2)
671 v0sccor(i,j)=v0ijsccor
677 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
678 c write (iout,*) 'ntortyp',ntortyp
680 cc maxinter is maximum interaction sites
684 read (isccor,*,end=113,err=113)
685 & nterm_sccor(i,j),nlor_sccor(i,j)
689 do k=1,nterm_sccor(i,j)
690 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
692 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
695 do k=1,nlor_sccor(i,j)
696 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
697 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
698 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
699 &(1+vlor3sccor(k,i,j)**2)
701 v0sccor(i,j)=v0ijsccor
709 write (iout,'(/a/)') 'Torsional constants:'
712 write (iout,*) 'ityp',i,' jtyp',j
713 write (iout,*) 'Fourier constants'
714 do k=1,nterm_sccor(i,j)
715 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
717 write (iout,*) 'Lorenz constants'
718 do k=1,nlor_sccor(i,j)
719 write (iout,'(3(1pe15.5))')
720 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
727 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
728 C interaction energy of the Gly, Ala, and Pro prototypes.
732 write (iout,*) "Coefficients of the cumulants"
734 read (ifourier,*) nloctyp
736 read (ifourier,*,end=115,err=115)
737 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
739 write (iout,*) 'Type',i
740 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
782 c Ctilde(1,1,i)=0.0d0
783 c Ctilde(1,2,i)=0.0d0
784 c Ctilde(2,1,i)=0.0d0
785 c Ctilde(2,2,i)=0.0d0
807 c Dtilde(1,1,i)=0.0d0
808 c Dtilde(1,2,i)=0.0d0
809 c Dtilde(2,1,i)=0.0d0
810 c Dtilde(2,2,i)=0.0d0
811 EE(1,1,i)= b(10)+b(11)
812 EE(2,2,i)=-b(10)+b(11)
813 EE(2,1,i)= b(12)-b(13)
814 EE(1,2,i)= b(12)+b(13)
815 EE(1,1,-i)= b(10)+b(11)
816 EE(2,2,-i)=-b(10)+b(11)
817 EE(2,1,-i)=-b(12)+b(13)
818 EE(1,2,-i)=-b(12)-b(13)
824 c ee(2,1,i)=ee(1,2,i)
828 write (iout,*) 'Type',i
830 write(iout,*) B1(1,i),B1(2,i)
832 write(iout,*) B2(1,i),B2(2,i)
835 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
839 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
843 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
849 C Read electrostatic-interaction parameters
853 write (iout,'(/a)') 'Electrostatic interaction constants:'
854 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
855 & 'IT','JT','APP','BPP','AEL6','AEL3'
857 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
858 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
859 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
860 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
865 app (i,j)=epp(i,j)*rri*rri
866 bpp (i,j)=-2.0D0*epp(i,j)*rri
867 ael6(i,j)=elpp6(i,j)*4.2D0**6
868 ael3(i,j)=elpp3(i,j)*4.2D0**3
869 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
870 & ael6(i,j),ael3(i,j)
874 C Read side-chain interaction parameters.
876 read (isidep,*,end=117,err=117) ipot,expon
877 if (ipot.lt.1 .or. ipot.gt.5) then
878 write (iout,'(2a)') 'Error while reading SC interaction',
879 & 'potential file - unknown potential type.'
881 call MPI_Finalize(Ierror)
887 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
888 & ', exponents are ',expon,2*expon
889 goto (10,20,30,30,40) ipot
890 C----------------------- LJ potential ---------------------------------
891 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
892 & (sigma0(i),i=1,ntyp)
894 write (iout,'(/a/)') 'Parameters of the LJ potential:'
895 write (iout,'(a/)') 'The epsilon array:'
896 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
897 write (iout,'(/a)') 'One-body parameters:'
898 write (iout,'(a,4x,a)') 'residue','sigma'
899 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
902 C----------------------- LJK potential --------------------------------
903 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
904 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
906 write (iout,'(/a/)') 'Parameters of the LJK potential:'
907 write (iout,'(a/)') 'The epsilon array:'
908 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
909 write (iout,'(/a)') 'One-body parameters:'
910 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
911 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
915 C---------------------- GB or BP potential -----------------------------
916 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
917 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
919 C For the GB potential convert sigma'**2 into chi'
922 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
926 write (iout,'(/a/)') 'Parameters of the BP potential:'
927 write (iout,'(a/)') 'The epsilon array:'
928 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
929 write (iout,'(/a)') 'One-body parameters:'
930 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
932 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
933 & chip(i),alp(i),i=1,ntyp)
936 C--------------------- GBV potential -----------------------------------
937 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
938 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
939 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
941 write (iout,'(/a/)') 'Parameters of the GBV potential:'
942 write (iout,'(a/)') 'The epsilon array:'
943 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
944 write (iout,'(/a)') 'One-body parameters:'
945 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
946 & 's||/s_|_^2',' chip ',' alph '
947 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
948 & sigii(i),chip(i),alp(i),i=1,ntyp)
952 C-----------------------------------------------------------------------
953 C Calculate the "working" parameters of SC interactions.
961 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
962 sigma(j,i)=sigma(i,j)
963 rs0(i,j)=dwa16*sigma(i,j)
967 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
968 & 'Working parameters of the SC interactions:',
969 & ' a ',' b ',' augm ',' sigma ',' r0 ',
974 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
983 sigeps=dsign(1.0D0,epsij)
985 aa(i,j)=epsij*rrij*rrij
986 bb(i,j)=-sigeps*epsij*rrij
994 ratsig1=sigt2sq/sigt1sq
995 ratsig2=1.0D0/ratsig1
996 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
997 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
998 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1002 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1003 sigmaii(i,j)=rsum_max
1004 sigmaii(j,i)=rsum_max
1006 c sigmaii(i,j)=r0(i,j)
1007 c sigmaii(j,i)=r0(i,j)
1009 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1010 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1011 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1012 augm(i,j)=epsij*r_augm**(2*expon)
1013 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1020 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1021 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1022 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1028 C Define the SC-p interaction constants (hard-coded; old style)
1031 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1033 c aad(i,1)=0.3D0*4.0D0**12
1034 C Following line for constants currently implemented
1035 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1036 aad(i,1)=1.5D0*4.0D0**12
1037 c aad(i,1)=0.17D0*5.6D0**12
1039 C "Soft" SC-p repulsion
1041 C Following line for constants currently implemented
1042 c aad(i,1)=0.3D0*4.0D0**6
1043 C "Hard" SC-p repulsion
1044 bad(i,1)=3.0D0*4.0D0**6
1045 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1054 C 8/9/01 Read the SC-p interaction constants from file
1057 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1060 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1061 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1062 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1063 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1067 write (iout,*) "Parameters of SC-p interactions:"
1069 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1070 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1075 C Define the constants of the disulfide bridge
1079 c Old arbitrary potential - commented out.
1084 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1085 c energy surface of diethyl disulfide.
1086 c A. Liwo and U. Kozlowska, 11/24/03
1103 write (iout,'(/a)') "Disulfide bridge parameters:"
1104 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1105 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1106 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1107 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1111 111 write (iout,*) "Error reading bending energy parameters."
1113 112 write (iout,*) "Error reading rotamer energy parameters."
1115 113 write (iout,*) "Error reading torsional energy parameters."
1117 114 write (iout,*) "Error reading double torsional energy parameters."
1120 & "Error reading cumulant (multibody energy) parameters."
1122 116 write (iout,*) "Error reading electrostatic energy parameters."
1124 117 write (iout,*) "Error reading side chain interaction parameters."
1126 118 write (iout,*) "Error reading SCp interaction parameters."
1128 119 write (iout,*) "Error reading SCCOR parameters"
1131 call MPI_Finalize(Ierror)
1138 subroutine getenv_loc(var, val)
1139 character(*) var, val
1142 character(2000) line
1145 open (196,file='env',status='old',readonly,shared)
1147 c write(*,*)'looking for ',var
1148 10 read(196,*,err=11,end=11)line
1149 iread=index(line,var)
1150 c write(*,*)iread,' ',var,' ',line
1151 if (iread.eq.0) go to 10
1152 c write(*,*)'---> ',line
1158 iread=iread+ilen(var)+1
1159 read (line(iread:),*,err=12,end=12) val
1160 c write(*,*)'OK: ',var,' = ',val
1166 #elif (defined CRAY)
1167 integer lennam,lenval,ierror
1169 c getenv using a POSIX call, useful on the T3D
1170 c Sept 1996, comment out error check on advice of H. Pritchard
1173 if(lennam.le.0) stop '--error calling getenv--'
1174 call pxfgetenv(var,lennam,val,lenval,ierror)
1175 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1177 call getenv(var,val)