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_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
319 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
320 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
321 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
324 write (2,*) "End reading THETA_PDB"
330 C Read the parameters of the probability distribution/energy expression
331 C of the side chains.
334 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
338 dsc_inv(i)=1.0D0/dsc(i)
349 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
350 & ((blower(k,l,1),l=1,k),k=1,3)
351 censc(1,1,-i)=censc(1,1,i)
352 censc(2,1,-i)=censc(2,1,i)
353 censc(3,1,-i)=-censc(3,1,i)
356 read (irotam,*,end=112,err=112) bsc(j,i)
357 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
358 & ((blower(k,l,j),l=1,k),k=1,3)
359 censc(1,j,-i)=censc(1,j,i)
360 censc(2,j,-i)=censc(2,j,i)
361 censc(3,j,-i)=-censc(3,j,i)
362 C BSC is amplitude of Gaussian
369 akl=akl+blower(k,m,j)*blower(l,m,j)
373 if (((k.eq.3).and.(l.ne.3))
374 & .or.((l.eq.3).and.(k.ne.3))) then
375 gaussc(k,l,j,-i)=-akl
376 gaussc(l,k,j,-i)=-akl
388 write (iout,'(/a)') 'Parameters of side-chain local geometry'
393 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
394 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
395 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
396 & 'log h',(bsc(j,i),j=1,nlobi)
397 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
398 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
400 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
401 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
404 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
405 write (iout,'(a,f10.4,4(16x,f10.4))')
406 & 'Center ',(bsc(j,i),j=1,nlobi)
407 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
416 C Read scrot parameters for potentials determined from all-atom AM1 calculations
417 C added by Urszula Kozlowska 07/11/2007
420 read (irotam,*,end=112,err=112)
422 read (irotam,*,end=112,err=112)
425 read(irotam,*,end=112,err=112) sc_parmin(j,i)
430 C Read the parameters of the probability distribution/energy expression
431 C of the side chains.
434 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
438 dsc_inv(i)=1.0D0/dsc(i)
449 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
450 & ((blower(k,l,1),l=1,k),k=1,3)
452 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
453 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
454 & ((blower(k,l,j),l=1,k),k=1,3)
461 akl=akl+blower(k,m,j)*blower(l,m,j)
476 C Read torsional parameters in old format
478 read (itorp,*,end=113,err=113) ntortyp,nterm_old
479 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
480 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
485 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
491 write (iout,'(/a/)') 'Torsional constants:'
494 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
495 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
501 C Read torsional parameters
503 read (itorp,*,end=113,err=113) ntortyp
504 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
507 itortyp(i)=-itortyp(-i)
509 c write (iout,*) 'ntortyp',ntortyp
511 do j=-ntortyp+1,ntortyp-1
512 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
514 nterm(-i,-j,iblock)=nterm(i,j,iblock)
515 nlor(-i,-j,iblock)=nlor(i,j,iblock)
518 do k=1,nterm(i,j,iblock)
519 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
521 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
522 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
523 v0ij=v0ij+si*v1(k,i,j,iblock)
525 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
526 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
527 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
529 do k=1,nlor(i,j,iblock)
530 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
531 & vlor2(k,i,j),vlor3(k,i,j)
532 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
535 v0(-i,-j,iblock)=v0ij
541 write (iout,'(/a/)') 'Torsional constants:'
544 write (iout,*) 'ityp',i,' jtyp',j
545 write (iout,*) 'Fourier constants'
546 do k=1,nterm(i,j,iblock)
547 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
550 write (iout,*) 'Lorenz constants'
551 do k=1,nlor(i,j,iblock)
552 write (iout,'(3(1pe15.5))')
553 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
559 C 6/23/01 Read parameters for double torsionals
563 do j=-ntortyp+1,ntortyp-1
564 do k=-ntortyp+1,ntortyp-1
565 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
566 write (iout,*) "OK onelett",
569 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
570 & .or. t3.ne.toronelet(k)) then
571 write (iout,*) "Error in double torsional parameter file",
574 call MPI_Finalize(Ierror)
576 stop "Error in double torsional parameter file"
578 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
579 & ntermd_2(i,j,k,iblock)
580 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
581 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
582 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
583 & ntermd_1(i,j,k,iblock))
584 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
585 & ntermd_1(i,j,k,iblock))
586 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
587 & ntermd_1(i,j,k,iblock))
588 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
589 & ntermd_1(i,j,k,iblock))
590 C Martix of D parameters for one dimesional foureir series
591 do l=1,ntermd_1(i,j,k,iblock)
592 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
593 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
594 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
595 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
596 c write(iout,*) "whcodze" ,
597 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
599 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
600 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
601 & v2s(m,l,i,j,k,iblock),
602 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
603 C Martix of D parameters for two dimesional fourier series
604 do l=1,ntermd_2(i,j,k,iblock)
606 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
607 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
608 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
609 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
618 write (iout,*) 'Constants for double torsionals'
621 do j=-ntortyp,ntortyp
622 do k=-ntortyp,ntortyp
623 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
624 & ' nsingle',ntermd_1(i,j,k,iblock),
625 & ' ndouble',ntermd_2(i,j,k,iblock)
627 write (iout,*) 'Single angles:'
628 do l=1,ntermd_1(i,j,k,iblock)
629 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
630 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
631 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
632 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
635 write (iout,*) 'Pairs of angles:'
636 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
637 do l=1,ntermd_2(i,j,k,iblock)
638 write (iout,'(i5,20f10.5)')
639 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
642 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
643 do l=1,ntermd_2(i,j,k,iblock)
644 write (iout,'(i5,20f10.5)')
645 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
646 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
655 C Read of Side-chain backbone correlation parameters
656 C Modified 11 May 2012 by Adasko
659 read (isccor,*,end=113,err=113) nsccortyp
661 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
663 read (isccor,*,end=113,err=113) (isccortyp(i),i=-ntyp,ntyp)
665 c write (iout,*) 'ntortyp',ntortyp
667 cc maxinter is maximum interaction sites
671 read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
675 do k=1,nterm_sccor(i,j)
676 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
678 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
681 do k=1,nlor_sccor(i,j)
682 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
683 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
684 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
685 &(1+vlor3sccor(k,i,j)**2)
687 v0sccor(i,j)=v0ijsccor
694 write (iout,'(/a/)') 'Torsional constants:'
697 write (iout,*) 'ityp',i,' jtyp',j
698 write (iout,*) 'Fourier constants'
699 do k=1,nterm_sccor(i,j)
700 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
702 write (iout,*) 'Lorenz constants'
703 do k=1,nlor_sccor(i,j)
704 write (iout,'(3(1pe15.5))')
705 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
712 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
713 C interaction energy of the Gly, Ala, and Pro prototypes.
717 write (iout,*) "Coefficients of the cumulants"
719 read (ifourier,*) nloctyp
721 read (ifourier,*,end=115,err=115)
722 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
724 write (iout,*) 'Type',i
725 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
767 c Ctilde(1,1,i)=0.0d0
768 c Ctilde(1,2,i)=0.0d0
769 c Ctilde(2,1,i)=0.0d0
770 c Ctilde(2,2,i)=0.0d0
792 c Dtilde(1,1,i)=0.0d0
793 c Dtilde(1,2,i)=0.0d0
794 c Dtilde(2,1,i)=0.0d0
795 c Dtilde(2,2,i)=0.0d0
796 EE(1,1,i)= b(10)+b(11)
797 EE(2,2,i)=-b(10)+b(11)
798 EE(2,1,i)= b(12)-b(13)
799 EE(1,2,i)= b(12)+b(13)
800 EE(1,1,-i)= b(10)+b(11)
801 EE(2,2,-i)=-b(10)+b(11)
802 EE(2,1,-i)=-b(12)+b(13)
803 EE(1,2,-i)=-b(12)-b(13)
809 c ee(2,1,i)=ee(1,2,i)
813 write (iout,*) 'Type',i
815 write(iout,*) B1(1,i),B1(2,i)
817 write(iout,*) B2(1,i),B2(2,i)
820 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
824 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
828 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
833 C Read electrostatic-interaction parameters
837 write (iout,'(/a)') 'Electrostatic interaction constants:'
838 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
839 & 'IT','JT','APP','BPP','AEL6','AEL3'
841 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
842 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
843 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
844 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
849 app (i,j)=epp(i,j)*rri*rri
850 bpp (i,j)=-2.0D0*epp(i,j)*rri
851 ael6(i,j)=elpp6(i,j)*4.2D0**6
852 ael3(i,j)=elpp3(i,j)*4.2D0**3
853 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
854 & ael6(i,j),ael3(i,j)
858 C Read side-chain interaction parameters.
860 read (isidep,*,end=117,err=117) ipot,expon
861 if (ipot.lt.1 .or. ipot.gt.5) then
862 write (iout,'(2a)') 'Error while reading SC interaction',
863 & 'potential file - unknown potential type.'
865 call MPI_Finalize(Ierror)
871 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
872 & ', exponents are ',expon,2*expon
873 goto (10,20,30,30,40) ipot
874 C----------------------- LJ potential ---------------------------------
875 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
876 & (sigma0(i),i=1,ntyp)
878 write (iout,'(/a/)') 'Parameters of the LJ potential:'
879 write (iout,'(a/)') 'The epsilon array:'
880 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
881 write (iout,'(/a)') 'One-body parameters:'
882 write (iout,'(a,4x,a)') 'residue','sigma'
883 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
886 C----------------------- LJK potential --------------------------------
887 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
888 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
890 write (iout,'(/a/)') 'Parameters of the LJK potential:'
891 write (iout,'(a/)') 'The epsilon array:'
892 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
893 write (iout,'(/a)') 'One-body parameters:'
894 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
895 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
899 C---------------------- GB or BP potential -----------------------------
900 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
901 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
903 C For the GB potential convert sigma'**2 into chi'
906 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
910 write (iout,'(/a/)') 'Parameters of the BP potential:'
911 write (iout,'(a/)') 'The epsilon array:'
912 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
913 write (iout,'(/a)') 'One-body parameters:'
914 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
916 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
917 & chip(i),alp(i),i=1,ntyp)
920 C--------------------- GBV potential -----------------------------------
921 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
922 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
923 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
925 write (iout,'(/a/)') 'Parameters of the GBV potential:'
926 write (iout,'(a/)') 'The epsilon array:'
927 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
928 write (iout,'(/a)') 'One-body parameters:'
929 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
930 & 's||/s_|_^2',' chip ',' alph '
931 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
932 & sigii(i),chip(i),alp(i),i=1,ntyp)
936 C-----------------------------------------------------------------------
937 C Calculate the "working" parameters of SC interactions.
945 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
946 sigma(j,i)=sigma(i,j)
947 rs0(i,j)=dwa16*sigma(i,j)
951 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
952 & 'Working parameters of the SC interactions:',
953 & ' a ',' b ',' augm ',' sigma ',' r0 ',
958 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
967 sigeps=dsign(1.0D0,epsij)
969 aa(i,j)=epsij*rrij*rrij
970 bb(i,j)=-sigeps*epsij*rrij
978 ratsig1=sigt2sq/sigt1sq
979 ratsig2=1.0D0/ratsig1
980 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
981 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
982 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
986 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
987 sigmaii(i,j)=rsum_max
988 sigmaii(j,i)=rsum_max
990 c sigmaii(i,j)=r0(i,j)
991 c sigmaii(j,i)=r0(i,j)
993 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
994 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
995 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
996 augm(i,j)=epsij*r_augm**(2*expon)
997 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1004 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1005 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1006 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1012 C Define the SC-p interaction constants (hard-coded; old style)
1015 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1017 c aad(i,1)=0.3D0*4.0D0**12
1018 C Following line for constants currently implemented
1019 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1020 aad(i,1)=1.5D0*4.0D0**12
1021 c aad(i,1)=0.17D0*5.6D0**12
1023 C "Soft" SC-p repulsion
1025 C Following line for constants currently implemented
1026 c aad(i,1)=0.3D0*4.0D0**6
1027 C "Hard" SC-p repulsion
1028 bad(i,1)=3.0D0*4.0D0**6
1029 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1038 C 8/9/01 Read the SC-p interaction constants from file
1041 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1044 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1045 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1046 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1047 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1051 write (iout,*) "Parameters of SC-p interactions:"
1053 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1054 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1059 C Define the constants of the disulfide bridge
1063 c Old arbitrary potential - commented out.
1068 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1069 c energy surface of diethyl disulfide.
1070 c A. Liwo and U. Kozlowska, 11/24/03
1087 write (iout,'(/a)') "Disulfide bridge parameters:"
1088 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1089 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1090 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1091 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1095 111 write (iout,*) "Error reading bending energy parameters."
1097 112 write (iout,*) "Error reading rotamer energy parameters."
1099 113 write (iout,*) "Error reading torsional energy parameters."
1101 114 write (iout,*) "Error reading double torsional energy parameters."
1104 & "Error reading cumulant (multibody energy) parameters."
1106 116 write (iout,*) "Error reading electrostatic energy parameters."
1108 117 write (iout,*) "Error reading side chain interaction parameters."
1110 118 write (iout,*) "Error reading SCp interaction parameters."
1112 119 write (iout,*) "Error reading SCCOR parameters"
1115 call MPI_Finalize(Ierror)
1122 subroutine getenv_loc(var, val)
1123 character(*) var, val
1126 character(2000) line
1129 open (196,file='env',status='old',readonly,shared)
1131 c write(*,*)'looking for ',var
1132 10 read(196,*,err=11,end=11)line
1133 iread=index(line,var)
1134 c write(*,*)iread,' ',var,' ',line
1135 if (iread.eq.0) go to 10
1136 c write(*,*)'---> ',line
1142 iread=iread+ilen(var)+1
1143 read (line(iread:),*,err=12,end=12) val
1144 c write(*,*)'OK: ',var,' = ',val
1150 #elif (defined CRAY)
1151 integer lennam,lenval,ierror
1153 c getenv using a POSIX call, useful on the T3D
1154 c Sept 1996, comment out error check on advice of H. Pritchard
1157 if(lennam.le.0) stop '--error calling getenv--'
1158 call pxfgetenv(var,lennam,val,lenval,ierror)
1159 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1161 call getenv(var,val)