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)
151 & 'Parameters of the virtual-bond valence angles:'
152 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
153 & ' ATHETA0 ',' A1 ',' A2 ',
156 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
157 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
159 write (iout,'(/a/9x,5a/79(1h-))')
160 & 'Parameters of the expression for sigma(theta_c):',
161 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
162 & ' ALPH3 ',' SIGMA0C '
164 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
165 & (polthet(j,i),j=0,3),sigc0(i)
167 write (iout,'(/a/9x,5a/79(1h-))')
168 & 'Parameters of the second gaussian:',
169 & ' THETA0 ',' SIGMA0 ',' G1 ',
172 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
173 & sig0(i),(gthet(j,i),j=1,3)
177 & 'Parameters of the virtual-bond valence angles:'
178 write (iout,'(/a/9x,5a/79(1h-))')
179 & 'Coefficients of expansion',
180 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
181 & ' b1*10^1 ',' b2*10^1 '
183 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
184 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
185 & (10*bthet(j,i,1,1),j=1,2)
187 write (iout,'(/a/9x,5a/79(1h-))')
188 & 'Parameters of the expression for sigma(theta_c):',
189 & ' alpha0 ',' alph1 ',' alph2 ',
190 & ' alhp3 ',' sigma0c '
192 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
193 & (polthet(j,i),j=0,3),sigc0(i)
195 write (iout,'(/a/9x,5a/79(1h-))')
196 & 'Parameters of the second gaussian:',
197 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
200 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
201 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
207 C Read the parameters of Utheta determined from ab initio surfaces
208 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
210 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
211 & ntheterm3,nsingle,ndouble
212 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
213 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
219 aathet(l,i,j,k)=0.0d0
223 bbthet(m,l,i,j,k)=0.0d0
224 ccthet(m,l,i,j,k)=0.0d0
225 ddthet(m,l,i,j,k)=0.0d0
226 eethet(m,l,i,j,k)=0.0d0
232 ffthet(mm,m,l,i,j,k)=0.0d0
233 ggthet(mm,m,l,i,j,k)=0.0d0
243 read (ithep,'(3a)',end=111,err=111) res1,res2,res3
244 read (ithep,*,end=111,err=111) aa0thet(i,j,k)
245 read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
246 read (ithep,*,end=111,err=111)
247 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
248 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
249 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
250 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
251 read (ithep,*,end=111,err=111)
252 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
253 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
254 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
259 C For dummy ends assign glycine-type coefficients of theta-only terms; the
260 C coefficients of theta-and-gamma-dependent terms are zero.
265 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
266 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
268 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
269 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
272 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
274 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
277 C Control printout of the coefficients of virtual-bond-angle potentials
280 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
284 write (iout,'(//4a)')
285 & 'Type ',onelett(i),onelett(j),onelett(k)
286 write (iout,'(//a,10x,a)') " l","a[l]"
287 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
288 write (iout,'(i2,1pe15.5)')
289 & (l,aathet(l,i,j,k),l=1,ntheterm)
291 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
292 & "b",l,"c",l,"d",l,"e",l
294 write (iout,'(i2,4(1pe15.5))') m,
295 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
296 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
300 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
301 & "f+",l,"f-",l,"g+",l,"g-",l
304 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
305 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
306 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
315 write (2,*) "Start reading THETA_PDB"
317 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
318 & (bthet(j,i,1,1),j=1,2)
319 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
320 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
321 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
325 athet(1,i,1,-1)=athet(1,i,1,1)
326 athet(2,i,1,-1)=athet(2,i,1,1)
327 bthet(1,i,1,-1)=-bthet(1,i,1,1)
328 bthet(2,i,1,-1)=-bthet(2,i,1,1)
329 athet(1,i,-1,1)=-athet(1,i,1,1)
330 athet(2,i,-1,1)=-athet(2,i,1,1)
331 bthet(1,i,-1,1)=bthet(1,i,1,1)
332 bthet(2,i,-1,1)=bthet(2,i,1,1)
336 athet(1,i,-1,-1)=athet(1,-i,1,1)
337 athet(2,i,-1,-1)=-athet(2,-i,1,1)
338 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
339 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
340 athet(1,i,-1,1)=athet(1,-i,1,1)
341 athet(2,i,-1,1)=-athet(2,-i,1,1)
342 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
343 bthet(2,i,-1,1)=bthet(2,-i,1,1)
344 athet(1,i,1,-1)=-athet(1,-i,1,1)
345 athet(2,i,1,-1)=athet(2,-i,1,1)
346 bthet(1,i,1,-1)=bthet(1,-i,1,1)
347 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
352 polthet(j,i)=polthet(j,-i)
355 gthet(j,i)=gthet(j,-i)
358 write (2,*) "End reading THETA_PDB"
364 C Read the parameters of the probability distribution/energy expression
365 C of the side chains.
368 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
372 dsc_inv(i)=1.0D0/dsc(i)
383 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
384 & ((blower(k,l,1),l=1,k),k=1,3)
385 censc(1,1,-i)=censc(1,1,i)
386 censc(2,1,-i)=censc(2,1,i)
387 censc(3,1,-i)=-censc(3,1,i)
390 read (irotam,*,end=112,err=112) bsc(j,i)
391 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
392 & ((blower(k,l,j),l=1,k),k=1,3)
393 censc(1,j,-i)=censc(1,j,i)
394 censc(2,j,-i)=censc(2,j,i)
395 censc(3,j,-i)=-censc(3,j,i)
396 C BSC is amplitude of Gaussian
403 akl=akl+blower(k,m,j)*blower(l,m,j)
407 if (((k.eq.3).and.(l.ne.3))
408 & .or.((l.eq.3).and.(k.ne.3))) then
409 gaussc(k,l,j,-i)=-akl
410 gaussc(l,k,j,-i)=-akl
422 write (iout,'(/a)') 'Parameters of side-chain local geometry'
427 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
428 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
429 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
430 & 'log h',(bsc(j,i),j=1,nlobi)
431 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
432 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
434 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
435 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
438 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
439 write (iout,'(a,f10.4,4(16x,f10.4))')
440 & 'Center ',(bsc(j,i),j=1,nlobi)
441 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
450 C Read scrot parameters for potentials determined from all-atom AM1 calculations
451 C added by Urszula Kozlowska 07/11/2007
454 read (irotam,*,end=112,err=112)
456 read (irotam,*,end=112,err=112)
459 read(irotam,*,end=112,err=112) sc_parmin(j,i)
464 C Read the parameters of the probability distribution/energy expression
465 C of the side chains.
468 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
472 dsc_inv(i)=1.0D0/dsc(i)
483 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
484 & ((blower(k,l,1),l=1,k),k=1,3)
486 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
487 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
488 & ((blower(k,l,j),l=1,k),k=1,3)
495 akl=akl+blower(k,m,j)*blower(l,m,j)
510 C Read torsional parameters in old format
512 read (itorp,*,end=113,err=113) ntortyp,nterm_old
513 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
514 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
519 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
525 write (iout,'(/a/)') 'Torsional constants:'
528 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
529 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
535 C Read torsional parameters
537 read (itorp,*,end=113,err=113) ntortyp
538 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
541 itortyp(i)=-itortyp(-i)
543 c write (iout,*) 'ntortyp',ntortyp
545 do j=-ntortyp+1,ntortyp-1
546 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
548 nterm(-i,-j,iblock)=nterm(i,j,iblock)
549 nlor(-i,-j,iblock)=nlor(i,j,iblock)
552 do k=1,nterm(i,j,iblock)
553 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
555 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
556 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
557 v0ij=v0ij+si*v1(k,i,j,iblock)
559 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
560 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
561 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
563 do k=1,nlor(i,j,iblock)
564 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
565 & vlor2(k,i,j),vlor3(k,i,j)
566 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
569 v0(-i,-j,iblock)=v0ij
575 write (iout,'(/a/)') 'Torsional constants:'
578 write (iout,*) 'ityp',i,' jtyp',j
579 write (iout,*) 'Fourier constants'
580 do k=1,nterm(i,j,iblock)
581 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
584 write (iout,*) 'Lorenz constants'
585 do k=1,nlor(i,j,iblock)
586 write (iout,'(3(1pe15.5))')
587 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
593 C 6/23/01 Read parameters for double torsionals
597 do j=-ntortyp+1,ntortyp-1
598 do k=-ntortyp+1,ntortyp-1
599 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
600 c write (iout,*) "OK onelett",
603 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
604 & .or. t3.ne.toronelet(k)) then
605 write (iout,*) "Error in double torsional parameter file",
608 call MPI_Finalize(Ierror)
610 stop "Error in double torsional parameter file"
612 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
613 & ntermd_2(i,j,k,iblock)
614 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
615 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
616 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
617 & ntermd_1(i,j,k,iblock))
618 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
619 & ntermd_1(i,j,k,iblock))
620 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
621 & ntermd_1(i,j,k,iblock))
622 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
623 & ntermd_1(i,j,k,iblock))
624 C Martix of D parameters for one dimesional foureir series
625 do l=1,ntermd_1(i,j,k,iblock)
626 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
627 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
628 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
629 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
630 c write(iout,*) "whcodze" ,
631 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
633 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
634 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
635 & v2s(m,l,i,j,k,iblock),
636 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
637 C Martix of D parameters for two dimesional fourier series
638 do l=1,ntermd_2(i,j,k,iblock)
640 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
641 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
642 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
643 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
652 write (iout,*) 'Constants for double torsionals'
655 do j=-ntortyp,ntortyp
656 do k=-ntortyp,ntortyp
657 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
658 & ' nsingle',ntermd_1(i,j,k,iblock),
659 & ' ndouble',ntermd_2(i,j,k,iblock)
661 write (iout,*) 'Single angles:'
662 do l=1,ntermd_1(i,j,k,iblock)
663 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
664 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
665 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
666 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
669 write (iout,*) 'Pairs of angles:'
670 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
671 do l=1,ntermd_2(i,j,k,iblock)
672 write (iout,'(i5,20f10.5)')
673 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
676 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
677 do l=1,ntermd_2(i,j,k,iblock)
678 write (iout,'(i5,20f10.5)')
679 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
680 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
689 C Read of Side-chain backbone correlation parameters
690 C Modified 11 May 2012 by Adasko
693 read (isccor,*,end=113,err=113) nsccortyp
695 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
697 isccortyp(i)=-isccortyp(-i)
699 iscprol=isccortyp(20)
700 c write (iout,*) 'ntortyp',ntortyp
702 cc maxinter is maximum interaction sites
706 read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
712 nterm_sccor(-i,j)=nterm_sccor(i,j)
713 nterm_sccor(-i,-j)=nterm_sccor(i,j)
714 nterm_sccor(i,-j)=nterm_sccor(i,j)
715 do k=1,nterm_sccor(i,j)
716 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
718 if (j.eq.iscprol) then
719 if (i.eq.isccortyp(10)) then
720 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
721 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
723 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
724 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
725 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
726 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
727 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
728 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
729 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
730 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
733 if (i.eq.isccortyp(10)) then
734 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
735 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
737 if (j.eq.isccortyp(10)) then
738 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
739 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
741 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
742 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
743 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
744 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
745 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
746 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
750 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
751 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
752 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
753 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
756 do k=1,nlor_sccor(i,j)
757 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
758 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
759 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
760 &(1+vlor3sccor(k,i,j)**2)
762 v0sccor(l,i,j)=v0ijsccor
763 v0sccor(l,-i,j)=v0ijsccor1
764 v0sccor(l,i,-j)=v0ijsccor2
765 v0sccor(l,-i,-j)=v0ijsccor3
771 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
772 c write (iout,*) 'ntortyp',ntortyp
774 cc maxinter is maximum interaction sites
778 read (isccor,*,end=113,err=113)
779 & nterm_sccor(i,j),nlor_sccor(i,j)
783 do k=1,nterm_sccor(i,j)
784 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
786 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
789 do k=1,nlor_sccor(i,j)
790 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
791 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
792 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
793 &(1+vlor3sccor(k,i,j)**2)
795 v0sccor(i,j)=v0ijsccor
803 write (iout,'(/a/)') 'Torsional constants:'
806 write (iout,*) 'ityp',i,' jtyp',j
807 write (iout,*) 'Fourier constants'
808 do k=1,nterm_sccor(i,j)
809 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
811 write (iout,*) 'Lorenz constants'
812 do k=1,nlor_sccor(i,j)
813 write (iout,'(3(1pe15.5))')
814 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
821 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
822 C interaction energy of the Gly, Ala, and Pro prototypes.
826 write (iout,*) "Coefficients of the cumulants"
828 read (ifourier,*) nloctyp
830 read (ifourier,*,end=115,err=115)
831 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
833 write (iout,*) 'Type',i
834 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
876 c Ctilde(1,1,i)=0.0d0
877 c Ctilde(1,2,i)=0.0d0
878 c Ctilde(2,1,i)=0.0d0
879 c Ctilde(2,2,i)=0.0d0
901 c Dtilde(1,1,i)=0.0d0
902 c Dtilde(1,2,i)=0.0d0
903 c Dtilde(2,1,i)=0.0d0
904 c Dtilde(2,2,i)=0.0d0
905 EE(1,1,i)= b(10)+b(11)
906 EE(2,2,i)=-b(10)+b(11)
907 EE(2,1,i)= b(12)-b(13)
908 EE(1,2,i)= b(12)+b(13)
909 EE(1,1,-i)= b(10)+b(11)
910 EE(2,2,-i)=-b(10)+b(11)
911 EE(2,1,-i)=-b(12)+b(13)
912 EE(1,2,-i)=-b(12)-b(13)
918 c ee(2,1,i)=ee(1,2,i)
922 write (iout,*) 'Type',i
924 write(iout,*) B1(1,i),B1(2,i)
926 write(iout,*) B2(1,i),B2(2,i)
929 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
933 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
937 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
942 C Read electrostatic-interaction parameters
946 write (iout,'(/a)') 'Electrostatic interaction constants:'
947 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
948 & 'IT','JT','APP','BPP','AEL6','AEL3'
950 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
951 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
952 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
953 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
958 app (i,j)=epp(i,j)*rri*rri
959 bpp (i,j)=-2.0D0*epp(i,j)*rri
960 ael6(i,j)=elpp6(i,j)*4.2D0**6
961 ael3(i,j)=elpp3(i,j)*4.2D0**3
962 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
963 & ael6(i,j),ael3(i,j)
967 C Read side-chain interaction parameters.
969 read (isidep,*,end=117,err=117) ipot,expon
970 if (ipot.lt.1 .or. ipot.gt.5) then
971 write (iout,'(2a)') 'Error while reading SC interaction',
972 & 'potential file - unknown potential type.'
974 call MPI_Finalize(Ierror)
980 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
981 & ', exponents are ',expon,2*expon
982 goto (10,20,30,30,40) ipot
983 C----------------------- LJ potential ---------------------------------
984 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
985 & (sigma0(i),i=1,ntyp)
987 write (iout,'(/a/)') 'Parameters of the LJ potential:'
988 write (iout,'(a/)') 'The epsilon array:'
989 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
990 write (iout,'(/a)') 'One-body parameters:'
991 write (iout,'(a,4x,a)') 'residue','sigma'
992 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
995 C----------------------- LJK potential --------------------------------
996 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
997 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
999 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1000 write (iout,'(a/)') 'The epsilon array:'
1001 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1002 write (iout,'(/a)') 'One-body parameters:'
1003 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1004 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1008 C---------------------- GB or BP potential -----------------------------
1009 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1010 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1012 C For the GB potential convert sigma'**2 into chi'
1015 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1019 write (iout,'(/a/)') 'Parameters of the BP 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,4a)') 'residue',' sigma ','s||/s_|_^2',
1025 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1026 & chip(i),alp(i),i=1,ntyp)
1029 C--------------------- GBV potential -----------------------------------
1030 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1031 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1032 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1034 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1035 write (iout,'(a/)') 'The epsilon array:'
1036 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1037 write (iout,'(/a)') 'One-body parameters:'
1038 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1039 & 's||/s_|_^2',' chip ',' alph '
1040 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1041 & sigii(i),chip(i),alp(i),i=1,ntyp)
1045 C-----------------------------------------------------------------------
1046 C Calculate the "working" parameters of SC interactions.
1054 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1055 sigma(j,i)=sigma(i,j)
1056 rs0(i,j)=dwa16*sigma(i,j)
1060 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1061 & 'Working parameters of the SC interactions:',
1062 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1067 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1076 sigeps=dsign(1.0D0,epsij)
1078 aa(i,j)=epsij*rrij*rrij
1079 bb(i,j)=-sigeps*epsij*rrij
1083 sigt1sq=sigma0(i)**2
1084 sigt2sq=sigma0(j)**2
1087 ratsig1=sigt2sq/sigt1sq
1088 ratsig2=1.0D0/ratsig1
1089 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1090 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1091 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1095 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1096 sigmaii(i,j)=rsum_max
1097 sigmaii(j,i)=rsum_max
1099 c sigmaii(i,j)=r0(i,j)
1100 c sigmaii(j,i)=r0(i,j)
1102 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1103 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1104 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1105 augm(i,j)=epsij*r_augm**(2*expon)
1106 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1113 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1114 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1115 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1121 C Define the SC-p interaction constants (hard-coded; old style)
1124 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1126 c aad(i,1)=0.3D0*4.0D0**12
1127 C Following line for constants currently implemented
1128 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1129 aad(i,1)=1.5D0*4.0D0**12
1130 c aad(i,1)=0.17D0*5.6D0**12
1132 C "Soft" SC-p repulsion
1134 C Following line for constants currently implemented
1135 c aad(i,1)=0.3D0*4.0D0**6
1136 C "Hard" SC-p repulsion
1137 bad(i,1)=3.0D0*4.0D0**6
1138 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1147 C 8/9/01 Read the SC-p interaction constants from file
1150 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1153 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1154 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1155 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1156 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1160 write (iout,*) "Parameters of SC-p interactions:"
1162 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1163 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1168 C Define the constants of the disulfide bridge
1172 c Old arbitrary potential - commented out.
1177 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1178 c energy surface of diethyl disulfide.
1179 c A. Liwo and U. Kozlowska, 11/24/03
1196 write (iout,'(/a)') "Disulfide bridge parameters:"
1197 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1198 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1199 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1200 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1204 111 write (iout,*) "Error reading bending energy parameters."
1206 112 write (iout,*) "Error reading rotamer energy parameters."
1208 113 write (iout,*) "Error reading torsional energy parameters."
1210 114 write (iout,*) "Error reading double torsional energy parameters."
1213 & "Error reading cumulant (multibody energy) parameters."
1215 116 write (iout,*) "Error reading electrostatic energy parameters."
1217 117 write (iout,*) "Error reading side chain interaction parameters."
1219 118 write (iout,*) "Error reading SCp interaction parameters."
1221 119 write (iout,*) "Error reading SCCOR parameters"
1224 call MPI_Finalize(Ierror)
1231 subroutine getenv_loc(var, val)
1232 character(*) var, val
1235 character(2000) line
1238 open (196,file='env',status='old',readonly,shared)
1240 c write(*,*)'looking for ',var
1241 10 read(196,*,err=11,end=11)line
1242 iread=index(line,var)
1243 c write(*,*)iread,' ',var,' ',line
1244 if (iread.eq.0) go to 10
1245 c write(*,*)'---> ',line
1251 iread=iread+ilen(var)+1
1252 read (line(iread:),*,err=12,end=12) val
1253 c write(*,*)'OK: ',var,' = ',val
1259 #elif (defined CRAY)
1260 integer lennam,lenval,ierror
1262 c getenv using a POSIX call, useful on the T3D
1263 c Sept 1996, comment out error check on advice of H. Pritchard
1266 if(lennam.le.0) stop '--error calling getenv--'
1267 call pxfgetenv(var,lennam,val,lenval,ierror)
1268 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1270 call getenv(var,val)