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"/
32 dimension blower(3,3,maxlob)
34 character*3 lancuch,ucase
36 C For printing parameters after they are read set the following in the UNRES
39 C setenv PRINT_PARM YES
41 C To print parameters in LaTeX format rather than as ASCII tables:
45 call getenv_loc("PRINT_PARM",lancuch)
46 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
47 call getenv_loc("LATEX",lancuch)
48 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
50 dwa16=2.0d0**(1.0d0/6.0d0)
52 C Assign virtual-bond length
57 c Read the virtual-bond parameters, masses, and moments of inertia
58 c and Stokes' radii of the peptide group and side chains
61 read (ibond,*) vbldp0,akp,mp,ip,pstok
64 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
69 dsc_inv(i)=1.0D0/dsc(i)
73 read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
75 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
76 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
81 dsc_inv(i)=1.0D0/dsc(i)
86 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
87 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
89 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
91 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
92 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
94 write (iout,'(13x,3f10.5)')
95 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
101 C Read the parameters of the probability distribution/energy expression
102 C of the virtual-bond valence angles theta
105 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
107 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
108 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
109 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
116 & 'Parameters of the virtual-bond valence angles:'
117 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
118 & ' ATHETA0 ',' A1 ',' A2 ',
121 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
122 & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
124 write (iout,'(/a/9x,5a/79(1h-))')
125 & 'Parameters of the expression for sigma(theta_c):',
126 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
127 & ' ALPH3 ',' SIGMA0C '
129 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
130 & (polthet(j,i),j=0,3),sigc0(i)
132 write (iout,'(/a/9x,5a/79(1h-))')
133 & 'Parameters of the second gaussian:',
134 & ' THETA0 ',' SIGMA0 ',' G1 ',
137 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
138 & sig0(i),(gthet(j,i),j=1,3)
142 & 'Parameters of the virtual-bond valence angles:'
143 write (iout,'(/a/9x,5a/79(1h-))')
144 & 'Coefficients of expansion',
145 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
146 & ' b1*10^1 ',' b2*10^1 '
148 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
149 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
151 write (iout,'(/a/9x,5a/79(1h-))')
152 & 'Parameters of the expression for sigma(theta_c):',
153 & ' alpha0 ',' alph1 ',' alph2 ',
154 & ' alhp3 ',' sigma0c '
156 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
157 & (polthet(j,i),j=0,3),sigc0(i)
159 write (iout,'(/a/9x,5a/79(1h-))')
160 & 'Parameters of the second gaussian:',
161 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
164 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
165 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
171 C Read the parameters of Utheta determined from ab initio surfaces
172 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
174 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
175 & ntheterm3,nsingle,ndouble
176 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
177 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
183 aathet(l,i,j,k)=0.0d0
187 bbthet(m,l,i,j,k)=0.0d0
188 ccthet(m,l,i,j,k)=0.0d0
189 ddthet(m,l,i,j,k)=0.0d0
190 eethet(m,l,i,j,k)=0.0d0
196 ffthet(mm,m,l,i,j,k)=0.0d0
197 ggthet(mm,m,l,i,j,k)=0.0d0
207 read (ithep,'(3a)',end=111,err=111) res1,res2,res3
208 read (ithep,*,end=111,err=111) aa0thet(i,j,k)
209 read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
210 read (ithep,*,end=111,err=111)
211 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
212 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
213 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
214 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
215 read (ithep,*,end=111,err=111)
216 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
217 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
218 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
223 C For dummy ends assign glycine-type coefficients of theta-only terms; the
224 C coefficients of theta-and-gamma-dependent terms are zero.
229 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
230 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
232 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
233 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
236 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
238 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
241 C Control printout of the coefficients of virtual-bond-angle potentials
244 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
248 write (iout,'(//4a)')
249 & 'Type ',onelett(i),onelett(j),onelett(k)
250 write (iout,'(//a,10x,a)') " l","a[l]"
251 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
252 write (iout,'(i2,1pe15.5)')
253 & (l,aathet(l,i,j,k),l=1,ntheterm)
255 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
256 & "b",l,"c",l,"d",l,"e",l
258 write (iout,'(i2,4(1pe15.5))') m,
259 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
260 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
264 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
265 & "f+",l,"f-",l,"g+",l,"g-",l
268 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
269 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
270 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
279 write (2,*) "Start reading THETA_PDB"
281 read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
283 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
284 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
285 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
288 write (2,*) "End reading THETA_PDB"
294 C Read the parameters of the probability distribution/energy expression
295 C of the side chains.
298 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
302 dsc_inv(i)=1.0D0/dsc(i)
313 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
314 & ((blower(k,l,1),l=1,k),k=1,3)
315 censc(1,1,-i)=censc(1,1,i)
316 censc(2,1,-i)=censc(2,1,i)
317 censc(3,1,-i)=-censc(3,1,i)
320 read (irotam,*,end=112,err=112) bsc(j,i)
321 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
322 & ((blower(k,l,j),l=1,k),k=1,3)
323 censc(1,j,-i)=censc(1,j,i)
324 censc(2,j,-i)=censc(2,j,i)
325 censc(3,j,-i)=-censc(3,j,i)
326 C BSC is amplitude of Gaussian
333 akl=akl+blower(k,m,j)*blower(l,m,j)
337 if (((k.eq.3).and.(l.ne.3))
338 & .or.((l.eq.3).and.(k.ne.3))) then
339 gaussc(k,l,j,-i)=-akl
340 gaussc(l,k,j,-i)=-akl
352 write (iout,'(/a)') 'Parameters of side-chain local geometry'
357 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
358 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
359 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
360 & 'log h',(bsc(j,i),j=1,nlobi)
361 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
362 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
364 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
365 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
368 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
369 write (iout,'(a,f10.4,4(16x,f10.4))')
370 & 'Center ',(bsc(j,i),j=1,nlobi)
371 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
380 C Read scrot parameters for potentials determined from all-atom AM1 calculations
381 C added by Urszula Kozlowska 07/11/2007
384 read (irotam,*,end=112,err=112)
386 read (irotam,*,end=112,err=112)
389 read(irotam,*,end=112,err=112) sc_parmin(j,i)
394 C Read the parameters of the probability distribution/energy expression
395 C of the side chains.
398 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
402 dsc_inv(i)=1.0D0/dsc(i)
413 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
414 & ((blower(k,l,1),l=1,k),k=1,3)
416 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
417 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
418 & ((blower(k,l,j),l=1,k),k=1,3)
425 akl=akl+blower(k,m,j)*blower(l,m,j)
440 C Read torsional parameters in old format
442 read (itorp,*,end=113,err=113) ntortyp,nterm_old
443 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
444 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
449 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
455 write (iout,'(/a/)') 'Torsional constants:'
458 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
459 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
465 C Read torsional parameters
467 read (itorp,*,end=113,err=113) ntortyp
468 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
471 itortyp(i)=-itortyp(-i)
473 c write (iout,*) 'ntortyp',ntortyp
475 do j=-ntortyp,ntortyp
476 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
478 nterm(-i,-j,iblock)=nterm(i,j,iblock)
479 nlor(-i,-j,iblock)=nlor(i,j,iblock)
482 do k=1,nterm(i,j,iblock)
483 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
485 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
486 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
487 v0ij=v0ij+si*v1(k,i,j)
490 do k=1,nlor(i,j,iblock)
491 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
492 & vlor2(k,i,j),vlor3(k,i,j)
493 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
496 v0(-i,-j,iblock)=v0ij
502 write (iout,'(/a/)') 'Torsional constants:'
505 write (iout,*) 'ityp',i,' jtyp',j
506 write (iout,*) 'Fourier constants'
508 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
510 write (iout,*) 'Lorenz constants'
512 write (iout,'(3(1pe15.5))')
513 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
519 C 6/23/01 Read parameters for double torsionals
524 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
525 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
526 & .or. t3.ne.onelett(k)) then
527 write (iout,*) "Error in double torsional parameter file",
530 call MPI_Finalize(Ierror)
532 stop "Error in double torsional parameter file"
534 read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
536 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
538 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
540 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
542 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
544 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
545 & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
546 & m=1,l-1),l=1,ntermd_2(i,j,k))
552 write (iout,*) 'Constants for double torsionals'
556 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
557 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
559 write (iout,*) 'Single angles:'
560 do l=1,ntermd_1(i,j,k)
561 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
562 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
563 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
566 write (iout,*) 'Pairs of angles:'
567 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
568 do l=1,ntermd_2(i,j,k)
569 write (iout,'(i5,20f10.5)')
570 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
573 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
574 do l=1,ntermd_2(i,j,k)
575 write (iout,'(i5,20f10.5)')
576 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
584 C Read of Side-chain backbone correlation parameters
585 C Modified 11 May 2012 by Adasko
588 read (isccor,*,end=113,err=113) nsccortyp
590 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
592 read (isccor,*,end=113,err=113) (isccortyp(i),i=-ntyp,ntyp)
594 c write (iout,*) 'ntortyp',ntortyp
596 cc maxinter is maximum interaction sites
600 read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
604 do k=1,nterm_sccor(i,j)
605 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
607 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
610 do k=1,nlor_sccor(i,j)
611 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
612 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
613 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
614 &(1+vlor3sccor(k,i,j)**2)
616 v0sccor(i,j)=v0ijsccor
623 write (iout,'(/a/)') 'Torsional constants:'
626 write (iout,*) 'ityp',i,' jtyp',j
627 write (iout,*) 'Fourier constants'
628 do k=1,nterm_sccor(i,j)
629 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
631 write (iout,*) 'Lorenz constants'
632 do k=1,nlor_sccor(i,j)
633 write (iout,'(3(1pe15.5))')
634 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
641 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
642 C interaction energy of the Gly, Ala, and Pro prototypes.
646 write (iout,*) "Coefficients of the cumulants"
648 read (ifourier,*) nloctyp
650 read (ifourier,*,end=115,err=115)
651 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
653 write (iout,*) 'Type',i
654 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
696 c Ctilde(1,1,i)=0.0d0
697 c Ctilde(1,2,i)=0.0d0
698 c Ctilde(2,1,i)=0.0d0
699 c Ctilde(2,2,i)=0.0d0
721 c Dtilde(1,1,i)=0.0d0
722 c Dtilde(1,2,i)=0.0d0
723 c Dtilde(2,1,i)=0.0d0
724 c Dtilde(2,2,i)=0.0d0
725 EE(1,1,i)= b(10)+b(11)
726 EE(2,2,i)=-b(10)+b(11)
727 EE(2,1,i)= b(12)-b(13)
728 EE(1,2,i)= b(12)+b(13)
729 EE(1,1,-i)= b(10)+b(11)
730 EE(2,2,-i)=-b(10)+b(11)
731 EE(2,1,-i)=-b(12)+b(13)
732 EE(1,2,-i)=-b(12)-b(13)
738 c ee(2,1,i)=ee(1,2,i)
742 write (iout,*) 'Type',i
744 write(iout,*) B1(1,i),B1(2,i)
746 write(iout,*) B2(1,i),B2(2,i)
749 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
753 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
757 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
762 C Read electrostatic-interaction parameters
766 write (iout,'(/a)') 'Electrostatic interaction constants:'
767 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
768 & 'IT','JT','APP','BPP','AEL6','AEL3'
770 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
771 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
772 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
773 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
778 app (i,j)=epp(i,j)*rri*rri
779 bpp (i,j)=-2.0D0*epp(i,j)*rri
780 ael6(i,j)=elpp6(i,j)*4.2D0**6
781 ael3(i,j)=elpp3(i,j)*4.2D0**3
782 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
783 & ael6(i,j),ael3(i,j)
787 C Read side-chain interaction parameters.
789 read (isidep,*,end=117,err=117) ipot,expon
790 if (ipot.lt.1 .or. ipot.gt.5) then
791 write (iout,'(2a)') 'Error while reading SC interaction',
792 & 'potential file - unknown potential type.'
794 call MPI_Finalize(Ierror)
800 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
801 & ', exponents are ',expon,2*expon
802 goto (10,20,30,30,40) ipot
803 C----------------------- LJ potential ---------------------------------
804 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
805 & (sigma0(i),i=1,ntyp)
807 write (iout,'(/a/)') 'Parameters of the LJ potential:'
808 write (iout,'(a/)') 'The epsilon array:'
809 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
810 write (iout,'(/a)') 'One-body parameters:'
811 write (iout,'(a,4x,a)') 'residue','sigma'
812 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
815 C----------------------- LJK potential --------------------------------
816 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
817 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
819 write (iout,'(/a/)') 'Parameters of the LJK potential:'
820 write (iout,'(a/)') 'The epsilon array:'
821 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
822 write (iout,'(/a)') 'One-body parameters:'
823 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
824 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
828 C---------------------- GB or BP potential -----------------------------
829 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
830 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
832 C For the GB potential convert sigma'**2 into chi'
835 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
839 write (iout,'(/a/)') 'Parameters of the BP potential:'
840 write (iout,'(a/)') 'The epsilon array:'
841 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
842 write (iout,'(/a)') 'One-body parameters:'
843 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
845 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
846 & chip(i),alp(i),i=1,ntyp)
849 C--------------------- GBV potential -----------------------------------
850 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
851 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
852 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
854 write (iout,'(/a/)') 'Parameters of the GBV potential:'
855 write (iout,'(a/)') 'The epsilon array:'
856 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
857 write (iout,'(/a)') 'One-body parameters:'
858 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
859 & 's||/s_|_^2',' chip ',' alph '
860 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
861 & sigii(i),chip(i),alp(i),i=1,ntyp)
865 C-----------------------------------------------------------------------
866 C Calculate the "working" parameters of SC interactions.
874 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
875 sigma(j,i)=sigma(i,j)
876 rs0(i,j)=dwa16*sigma(i,j)
880 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
881 & 'Working parameters of the SC interactions:',
882 & ' a ',' b ',' augm ',' sigma ',' r0 ',
887 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
896 sigeps=dsign(1.0D0,epsij)
898 aa(i,j)=epsij*rrij*rrij
899 bb(i,j)=-sigeps*epsij*rrij
907 ratsig1=sigt2sq/sigt1sq
908 ratsig2=1.0D0/ratsig1
909 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
910 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
911 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
915 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
916 sigmaii(i,j)=rsum_max
917 sigmaii(j,i)=rsum_max
919 c sigmaii(i,j)=r0(i,j)
920 c sigmaii(j,i)=r0(i,j)
922 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
923 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
924 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
925 augm(i,j)=epsij*r_augm**(2*expon)
926 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
933 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
934 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
935 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
941 C Define the SC-p interaction constants (hard-coded; old style)
944 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
946 c aad(i,1)=0.3D0*4.0D0**12
947 C Following line for constants currently implemented
948 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
949 aad(i,1)=1.5D0*4.0D0**12
950 c aad(i,1)=0.17D0*5.6D0**12
952 C "Soft" SC-p repulsion
954 C Following line for constants currently implemented
955 c aad(i,1)=0.3D0*4.0D0**6
956 C "Hard" SC-p repulsion
957 bad(i,1)=3.0D0*4.0D0**6
958 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
967 C 8/9/01 Read the SC-p interaction constants from file
970 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
973 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
974 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
975 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
976 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
980 write (iout,*) "Parameters of SC-p interactions:"
982 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
983 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
988 C Define the constants of the disulfide bridge
992 c Old arbitrary potential - commented out.
997 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
998 c energy surface of diethyl disulfide.
999 c A. Liwo and U. Kozlowska, 11/24/03
1016 write (iout,'(/a)') "Disulfide bridge parameters:"
1017 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1018 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1019 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1020 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1024 111 write (iout,*) "Error reading bending energy parameters."
1026 112 write (iout,*) "Error reading rotamer energy parameters."
1028 113 write (iout,*) "Error reading torsional energy parameters."
1030 114 write (iout,*) "Error reading double torsional energy parameters."
1033 & "Error reading cumulant (multibody energy) parameters."
1035 116 write (iout,*) "Error reading electrostatic energy parameters."
1037 117 write (iout,*) "Error reading side chain interaction parameters."
1039 118 write (iout,*) "Error reading SCp interaction parameters."
1041 119 write (iout,*) "Error reading SCCOR parameters"
1044 call MPI_Finalize(Ierror)
1051 subroutine getenv_loc(var, val)
1052 character(*) var, val
1055 character(2000) line
1058 open (196,file='env',status='old',readonly,shared)
1060 c write(*,*)'looking for ',var
1061 10 read(196,*,err=11,end=11)line
1062 iread=index(line,var)
1063 c write(*,*)iread,' ',var,' ',line
1064 if (iread.eq.0) go to 10
1065 c write(*,*)'---> ',line
1071 iread=iread+ilen(var)+1
1072 read (line(iread:),*,err=12,end=12) val
1073 c write(*,*)'OK: ',var,' = ',val
1079 #elif (defined CRAY)
1080 integer lennam,lenval,ierror
1082 c getenv using a POSIX call, useful on the T3D
1083 c Sept 1996, comment out error check on advice of H. Pritchard
1086 if(lennam.le.0) stop '--error calling getenv--'
1087 call pxfgetenv(var,lennam,val,lenval,ierror)
1088 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1090 call getenv(var,val)