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)
216 ithetyp(i)=-ithetyp(-i)
219 do i=-maxthetyp,maxthetyp
220 do j=-maxthetyp,maxthetyp
221 do k=-maxthetyp,maxthetyp
222 aa0thet(i,j,k,iblock)=0.0d0
224 aathet(l,i,j,k,iblock)=0.0d0
228 bbthet(m,l,i,j,k,iblock)=0.0d0
229 ccthet(m,l,i,j,k,iblock)=0.0d0
230 ddthet(m,l,i,j,k,iblock)=0.0d0
231 eethet(m,l,i,j,k,iblock)=0.0d0
237 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
238 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
248 do j=-nthetyp,nthetyp
249 do k=-nthetyp,nthetyp
250 read (ithep,'(6a)',end=111,err=111) res1
251 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
252 read (ithep,*,end=111,err=111)
253 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
254 read (ithep,*,end=111,err=111)
255 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
256 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
257 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
258 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
260 read (ithep,*,end=111,err=111)
261 & (((ffthet(llll,lll,ll,i,j,k,iblock),
262 & ffthet(lll,llll,ll,i,j,k,iblock),
263 & ggthet(llll,lll,ll,i,j,k,iblock),
264 & ggthet(lll,llll,ll,i,j,k,iblock),
265 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
270 C For dummy ends assign glycine-type coefficients of theta-only terms; the
271 C coefficients of theta-and-gamma-dependent terms are zero.
276 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
277 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
279 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
280 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
283 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
285 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
291 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
292 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
294 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
295 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
298 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
300 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
304 C Substitution for D aminoacids from symmetry.
307 do j=-nthetyp,nthetyp
308 do k=-nthetyp,nthetyp
309 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
311 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
315 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
316 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
317 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
318 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
324 ffthet(llll,lll,ll,i,j,k,iblock)=
325 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
326 ffthet(lll,llll,ll,i,j,k,iblock)=
327 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
328 ggthet(llll,lll,ll,i,j,k,iblock)=
329 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
330 ggthet(lll,llll,ll,i,j,k,iblock)=
331 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
340 C Control printout of the coefficients of virtual-bond-angle potentials
343 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
347 write (iout,'(//4a)')
348 & 'Type ',onelett(i),onelett(j),onelett(k)
349 write (iout,'(//a,10x,a)') " l","a[l]"
350 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
351 write (iout,'(i2,1pe15.5)')
352 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
354 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
355 & "b",l,"c",l,"d",l,"e",l
357 write (iout,'(i2,4(1pe15.5))') m,
358 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
359 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
363 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
364 & "f+",l,"f-",l,"g+",l,"g-",l
367 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
368 & ffthet(n,m,l,i,j,k,iblock),
369 & ffthet(m,n,l,i,j,k,iblock),
370 & ggthet(n,m,l,i,j,k,iblock),
371 & ggthet(m,n,l,i,j,k,iblock)
384 C Read the parameters of the probability distribution/energy expression
385 C of the side chains.
388 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
392 dsc_inv(i)=1.0D0/dsc(i)
403 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
404 & ((blower(k,l,1),l=1,k),k=1,3)
405 censc(1,1,-i)=censc(1,1,i)
406 censc(2,1,-i)=censc(2,1,i)
407 censc(3,1,-i)=-censc(3,1,i)
409 read (irotam,*,end=112,err=112) bsc(j,i)
410 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
411 & ((blower(k,l,j),l=1,k),k=1,3)
412 censc(1,j,-i)=censc(1,j,i)
413 censc(2,j,-i)=censc(2,j,i)
414 censc(3,j,-i)=-censc(3,j,i)
415 C BSC is amplitude of Gaussian
422 akl=akl+blower(k,m,j)*blower(l,m,j)
426 if (((k.eq.3).and.(l.ne.3))
427 & .or.((l.eq.3).and.(k.ne.3))) then
428 gaussc(k,l,j,-i)=-akl
429 gaussc(l,k,j,-i)=-akl
441 write (iout,'(/a)') 'Parameters of side-chain local geometry'
446 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
447 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
448 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
449 & 'log h',(bsc(j,i),j=1,nlobi)
450 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
451 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
453 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
454 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
457 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
458 write (iout,'(a,f10.4,4(16x,f10.4))')
459 & 'Center ',(bsc(j,i),j=1,nlobi)
460 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
469 C Read scrot parameters for potentials determined from all-atom AM1 calculations
470 C added by Urszula Kozlowska 07/11/2007
473 CC TU JEST ZLE musibyc ntyp
474 read (irotam,*,end=112,err=112)
476 read (irotam,*,end=112,err=112)
479 read(irotam,*,end=112,err=112) sc_parmin(j,i)
488 C Read torsional parameters in old format
490 read (itorp,*,end=113,err=113) ntortyp,nterm_old
491 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
492 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
497 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
503 write (iout,'(/a/)') 'Torsional constants:'
506 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
507 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
513 C Read torsional parameters
515 read (itorp,*,end=113,err=113) ntortyp
516 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
519 itortyp(i)=-itortyp(-i)
521 write (iout,*) 'ntortyp',ntortyp
523 do j=-ntortyp+1,ntortyp-1
524 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
526 nterm(-i,-j,iblock)=nterm(i,j,iblock)
527 nlor(-i,-j,iblock)=nlor(i,j,iblock)
530 do k=1,nterm(i,j,iblock)
531 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
533 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
534 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
535 v0ij=v0ij+si*v1(k,i,j,iblock)
537 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
538 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
539 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
541 do k=1,nlor(i,j,iblock)
542 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
543 & vlor2(k,i,j),vlor3(k,i,j)
544 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
547 v0(-i,-j,iblock)=v0ij
553 write (iout,'(/a/)') 'Torsional constants:'
556 write (iout,*) 'ityp',i,' jtyp',j
557 write (iout,*) 'Fourier constants'
558 do k=1,nterm(i,j,iblock)
559 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
562 write (iout,*) 'Lorenz constants'
563 do k=1,nlor(i,j,iblock)
564 write (iout,'(3(1pe15.5))')
565 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
572 C 6/23/01 Read parameters for double torsionals
576 do j=-ntortyp+1,ntortyp-1
577 do k=-ntortyp+1,ntortyp-1
578 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
579 c write (iout,*) "OK onelett",
582 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
583 & .or. t3.ne.toronelet(k)) then
584 write (iout,*) "Error in double torsional parameter file",
587 call MPI_Finalize(Ierror)
589 stop "Error in double torsional parameter file"
591 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
592 & ntermd_2(i,j,k,iblock)
593 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
594 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
595 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
596 & ntermd_1(i,j,k,iblock))
597 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
598 & ntermd_1(i,j,k,iblock))
599 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
600 & ntermd_1(i,j,k,iblock))
601 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
602 & ntermd_1(i,j,k,iblock))
603 C Martix of D parameters for one dimesional foureir series
604 do l=1,ntermd_1(i,j,k,iblock)
605 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
606 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
607 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
608 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
609 c write(iout,*) "whcodze" ,
610 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
612 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
613 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
614 & v2s(m,l,i,j,k,iblock),
615 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
616 C Martix of D parameters for two dimesional fourier series
617 do l=1,ntermd_2(i,j,k,iblock)
619 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
620 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
621 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
622 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
631 write (iout,*) 'Constants for double torsionals'
634 do j=-ntortyp+1,ntortyp-1
635 do k=-ntortyp+1,ntortyp-1
636 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
637 & ' nsingle',ntermd_1(i,j,k,iblock),
638 & ' ndouble',ntermd_2(i,j,k,iblock)
640 write (iout,*) 'Single angles:'
641 do l=1,ntermd_1(i,j,k,iblock)
642 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
643 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
644 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
645 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
648 write (iout,*) 'Pairs of angles:'
649 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
650 do l=1,ntermd_2(i,j,k,iblock)
651 write (iout,'(i5,20f10.5)')
652 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
655 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
656 do l=1,ntermd_2(i,j,k,iblock)
657 write (iout,'(i5,20f10.5)')
658 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
659 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
668 C Read of Side-chain backbone correlation parameters
669 C Modified 11 May 2012 by Adasko
672 read (isccor,*,end=119,err=119) nsccortyp
674 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
676 isccortyp(i)=-isccortyp(-i)
678 iscprol=isccortyp(20)
679 c write (iout,*) 'ntortyp',ntortyp
681 cc maxinter is maximum interaction sites
685 read (isccor,*,end=119,err=119)
686 &nterm_sccor(i,j),nlor_sccor(i,j)
692 nterm_sccor(-i,j)=nterm_sccor(i,j)
693 nterm_sccor(-i,-j)=nterm_sccor(i,j)
694 nterm_sccor(i,-j)=nterm_sccor(i,j)
695 do k=1,nterm_sccor(i,j)
696 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
698 if (j.eq.iscprol) then
699 if (i.eq.isccortyp(10)) then
700 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
701 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
703 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
704 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
705 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
706 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
707 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
708 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
709 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
710 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
713 if (i.eq.isccortyp(10)) then
714 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
715 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
717 if (j.eq.isccortyp(10)) then
718 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
719 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
721 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
722 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
723 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
724 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
725 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
726 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
730 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
731 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
732 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
733 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
736 do k=1,nlor_sccor(i,j)
737 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
738 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
739 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
740 &(1+vlor3sccor(k,i,j)**2)
742 v0sccor(l,i,j)=v0ijsccor
743 v0sccor(l,-i,j)=v0ijsccor1
744 v0sccor(l,i,-j)=v0ijsccor2
745 v0sccor(l,-i,-j)=v0ijsccor3
751 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
752 c write (iout,*) 'ntortyp',ntortyp
754 cc maxinter is maximum interaction sites
758 read (isccor,*,end=113,err=113)
759 & nterm_sccor(i,j),nlor_sccor(i,j)
763 do k=1,nterm_sccor(i,j)
764 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
766 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
769 do k=1,nlor_sccor(i,j)
770 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
771 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
772 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
773 &(1+vlor3sccor(k,i,j)**2)
775 v0sccor(i,j,iblock)=v0ijsccor
783 write (iout,'(/a/)') 'Torsional constants:'
786 write (iout,*) 'ityp',i,' jtyp',j
787 write (iout,*) 'Fourier constants'
788 do k=1,nterm_sccor(i,j)
789 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
791 write (iout,*) 'Lorenz constants'
792 do k=1,nlor_sccor(i,j)
793 write (iout,'(3(1pe15.5))')
794 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
801 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
802 C interaction energy of the Gly, Ala, and Pro prototypes.
806 write (iout,*) "Coefficients of the cumulants"
808 read (ifourier,*) nloctyp
810 read (ifourier,*,end=115,err=115)
811 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
813 write (iout,*) 'Type',i
814 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
856 c Ctilde(1,1,i)=0.0d0
857 c Ctilde(1,2,i)=0.0d0
858 c Ctilde(2,1,i)=0.0d0
859 c Ctilde(2,2,i)=0.0d0
881 c Dtilde(1,1,i)=0.0d0
882 c Dtilde(1,2,i)=0.0d0
883 c Dtilde(2,1,i)=0.0d0
884 c Dtilde(2,2,i)=0.0d0
885 EE(1,1,i)= b(10)+b(11)
886 EE(2,2,i)=-b(10)+b(11)
887 EE(2,1,i)= b(12)-b(13)
888 EE(1,2,i)= b(12)+b(13)
889 EE(1,1,-i)= b(10)+b(11)
890 EE(2,2,-i)=-b(10)+b(11)
891 EE(2,1,-i)=-b(12)+b(13)
892 EE(1,2,-i)=-b(12)-b(13)
898 c ee(2,1,i)=ee(1,2,i)
902 write (iout,*) 'Type',i
904 write(iout,*) B1(1,i),B1(2,i)
906 write(iout,*) B2(1,i),B2(2,i)
909 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
913 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
917 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
923 C Read electrostatic-interaction parameters
927 write (iout,'(/a)') 'Electrostatic interaction constants:'
928 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
929 & 'IT','JT','APP','BPP','AEL6','AEL3'
931 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
932 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
933 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
934 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
939 app (i,j)=epp(i,j)*rri*rri
940 bpp (i,j)=-2.0D0*epp(i,j)*rri
941 ael6(i,j)=elpp6(i,j)*4.2D0**6
942 ael3(i,j)=elpp3(i,j)*4.2D0**3
943 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
944 & ael6(i,j),ael3(i,j)
948 C Read side-chain interaction parameters.
950 read (isidep,*,end=117,err=117) ipot,expon
951 if (ipot.lt.1 .or. ipot.gt.5) then
952 write (iout,'(2a)') 'Error while reading SC interaction',
953 & 'potential file - unknown potential type.'
955 call MPI_Finalize(Ierror)
961 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
962 & ', exponents are ',expon,2*expon
963 goto (10,20,30,30,40) ipot
964 C----------------------- LJ potential ---------------------------------
965 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
966 & (sigma0(i),i=1,ntyp)
968 write (iout,'(/a/)') 'Parameters of the LJ potential:'
969 write (iout,'(a/)') 'The epsilon array:'
970 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
971 write (iout,'(/a)') 'One-body parameters:'
972 write (iout,'(a,4x,a)') 'residue','sigma'
973 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
976 C----------------------- LJK potential --------------------------------
977 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
978 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
980 write (iout,'(/a/)') 'Parameters of the LJK potential:'
981 write (iout,'(a/)') 'The epsilon array:'
982 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
983 write (iout,'(/a)') 'One-body parameters:'
984 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
985 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
989 C---------------------- GB or BP potential -----------------------------
991 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
993 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
994 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
995 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
996 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
997 C For the GB potential convert sigma'**2 into chi'
1000 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1004 write (iout,'(/a/)') 'Parameters of the BP potential:'
1005 write (iout,'(a/)') 'The epsilon array:'
1006 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1007 write (iout,'(/a)') 'One-body parameters:'
1008 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1010 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1011 & chip(i),alp(i),i=1,ntyp)
1014 C--------------------- GBV potential -----------------------------------
1015 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1016 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1017 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1019 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1020 write (iout,'(a/)') 'The epsilon array:'
1021 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1022 write (iout,'(/a)') 'One-body parameters:'
1023 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1024 & 's||/s_|_^2',' chip ',' alph '
1025 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1026 & sigii(i),chip(i),alp(i),i=1,ntyp)
1030 C-----------------------------------------------------------------------
1031 C Calculate the "working" parameters of SC interactions.
1039 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1040 sigma(j,i)=sigma(i,j)
1041 rs0(i,j)=dwa16*sigma(i,j)
1045 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1046 & 'Working parameters of the SC interactions:',
1047 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1052 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1061 sigeps=dsign(1.0D0,epsij)
1063 aa(i,j)=epsij*rrij*rrij
1064 bb(i,j)=-sigeps*epsij*rrij
1068 sigt1sq=sigma0(i)**2
1069 sigt2sq=sigma0(j)**2
1072 ratsig1=sigt2sq/sigt1sq
1073 ratsig2=1.0D0/ratsig1
1074 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1075 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1076 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1080 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1081 sigmaii(i,j)=rsum_max
1082 sigmaii(j,i)=rsum_max
1084 c sigmaii(i,j)=r0(i,j)
1085 c sigmaii(j,i)=r0(i,j)
1087 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1088 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1089 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1090 augm(i,j)=epsij*r_augm**(2*expon)
1091 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1098 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1099 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1100 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1106 C Define the SC-p interaction constants (hard-coded; old style)
1109 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1111 c aad(i,1)=0.3D0*4.0D0**12
1112 C Following line for constants currently implemented
1113 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1114 aad(i,1)=1.5D0*4.0D0**12
1115 c aad(i,1)=0.17D0*5.6D0**12
1117 C "Soft" SC-p repulsion
1119 C Following line for constants currently implemented
1120 c aad(i,1)=0.3D0*4.0D0**6
1121 C "Hard" SC-p repulsion
1122 bad(i,1)=3.0D0*4.0D0**6
1123 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1132 C 8/9/01 Read the SC-p interaction constants from file
1135 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1138 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1139 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1140 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1141 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1145 write (iout,*) "Parameters of SC-p interactions:"
1147 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1148 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1153 C Define the constants of the disulfide bridge
1157 c Old arbitrary potential - commented out.
1162 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1163 c energy surface of diethyl disulfide.
1164 c A. Liwo and U. Kozlowska, 11/24/03
1181 write (iout,'(/a)') "Disulfide bridge parameters:"
1182 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1183 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1184 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1185 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1189 111 write (iout,*) "Error reading bending energy parameters."
1191 112 write (iout,*) "Error reading rotamer energy parameters."
1193 113 write (iout,*) "Error reading torsional energy parameters."
1195 114 write (iout,*) "Error reading double torsional energy parameters."
1198 & "Error reading cumulant (multibody energy) parameters."
1200 116 write (iout,*) "Error reading electrostatic energy parameters."
1202 117 write (iout,*) "Error reading side chain interaction parameters."
1204 118 write (iout,*) "Error reading SCp interaction parameters."
1206 119 write (iout,*) "Error reading SCCOR parameters"
1209 call MPI_Finalize(Ierror)
1216 subroutine getenv_loc(var, val)
1217 character(*) var, val
1220 character(2000) line
1223 open (196,file='env',status='old',readonly,shared)
1225 c write(*,*)'looking for ',var
1226 10 read(196,*,err=11,end=11)line
1227 iread=index(line,var)
1228 c write(*,*)iread,' ',var,' ',line
1229 if (iread.eq.0) go to 10
1230 c write(*,*)'---> ',line
1236 iread=iread+ilen(var)+1
1237 read (line(iread:),*,err=12,end=12) val
1238 c write(*,*)'OK: ',var,' = ',val
1244 #elif (defined CRAY)
1245 integer lennam,lenval,ierror
1247 c getenv using a POSIX call, useful on the T3D
1248 c Sept 1996, comment out error check on advice of H. Pritchard
1251 if(lennam.le.0) stop '--error calling getenv--'
1252 call pxfgetenv(var,lennam,val,lenval,ierror)
1253 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1255 call getenv(var,val)