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),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)
117 & 'Parameters of the virtual-bond valence angles:'
118 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
119 & ' ATHETA0 ',' A1 ',' A2 ',
122 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
123 & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
125 write (iout,'(/a/9x,5a/79(1h-))')
126 & 'Parameters of the expression for sigma(theta_c):',
127 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
128 & ' ALPH3 ',' SIGMA0C '
130 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
131 & (polthet(j,i),j=0,3),sigc0(i)
133 write (iout,'(/a/9x,5a/79(1h-))')
134 & 'Parameters of the second gaussian:',
135 & ' THETA0 ',' SIGMA0 ',' G1 ',
138 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
139 & sig0(i),(gthet(j,i),j=1,3)
143 & 'Parameters of the virtual-bond valence angles:'
144 write (iout,'(/a/9x,5a/79(1h-))')
145 & 'Coefficients of expansion',
146 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
147 & ' b1*10^1 ',' b2*10^1 '
149 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
150 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
152 write (iout,'(/a/9x,5a/79(1h-))')
153 & 'Parameters of the expression for sigma(theta_c):',
154 & ' alpha0 ',' alph1 ',' alph2 ',
155 & ' alhp3 ',' sigma0c '
157 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
158 & (polthet(j,i),j=0,3),sigc0(i)
160 write (iout,'(/a/9x,5a/79(1h-))')
161 & 'Parameters of the second gaussian:',
162 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
165 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
166 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
172 C Read the parameters of Utheta determined from ab initio surfaces
173 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
175 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
176 & ntheterm3,nsingle,ndouble
177 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
178 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
184 aathet(l,i,j,k)=0.0d0
188 bbthet(m,l,i,j,k)=0.0d0
189 ccthet(m,l,i,j,k)=0.0d0
190 ddthet(m,l,i,j,k)=0.0d0
191 eethet(m,l,i,j,k)=0.0d0
197 ffthet(mm,m,l,i,j,k)=0.0d0
198 ggthet(mm,m,l,i,j,k)=0.0d0
208 read (ithep,'(3a)',end=111,err=111) res1,res2,res3
209 read (ithep,*,end=111,err=111) aa0thet(i,j,k)
210 read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
211 read (ithep,*,end=111,err=111)
212 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
213 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
214 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
215 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
216 read (ithep,*,end=111,err=111)
217 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
218 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
219 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
224 C For dummy ends assign glycine-type coefficients of theta-only terms; the
225 C coefficients of theta-and-gamma-dependent terms are zero.
230 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
231 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
233 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
234 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
237 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
239 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
242 C Control printout of the coefficients of virtual-bond-angle potentials
245 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
249 write (iout,'(//4a)')
250 & 'Type ',onelett(i),onelett(j),onelett(k)
251 write (iout,'(//a,10x,a)') " l","a[l]"
252 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
253 write (iout,'(i2,1pe15.5)')
254 & (l,aathet(l,i,j,k),l=1,ntheterm)
256 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
257 & "b",l,"c",l,"d",l,"e",l
259 write (iout,'(i2,4(1pe15.5))') m,
260 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
261 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
265 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
266 & "f+",l,"f-",l,"g+",l,"g-",l
269 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
270 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
271 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
280 write (2,*) "Start reading THETA_PDB"
282 read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
284 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
285 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
286 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
288 write (2,*) "End reading THETA_PDB"
295 C Read the parameters of the probability distribution/energy expression
296 C of the side chains.
299 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
303 dsc_inv(i)=1.0D0/dsc(i)
314 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
315 & ((blower(k,l,1),l=1,k),k=1,3)
317 read (irotam,*,end=112,err=112) bsc(j,i)
318 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
319 & ((blower(k,l,j),l=1,k),k=1,3)
326 akl=akl+blower(k,m,j)*blower(l,m,j)
337 write (iout,'(/a)') 'Parameters of side-chain local geometry'
342 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
343 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
344 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
345 & 'log h',(bsc(j,i),j=1,nlobi)
346 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
347 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
349 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
350 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
353 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
354 write (iout,'(a,f10.4,4(16x,f10.4))')
355 & 'Center ',(bsc(j,i),j=1,nlobi)
356 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
365 C Read scrot parameters for potentials determined from all-atom AM1 calculations
366 C added by Urszula Kozlowska 07/11/2007
369 read (irotam,*,end=112,err=112)
371 read (irotam,*,end=112,err=112)
374 read(irotam,*,end=112,err=112) sc_parmin(j,i)
379 C Read the parameters of the probability distribution/energy expression
380 C of the side chains.
383 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
387 dsc_inv(i)=1.0D0/dsc(i)
398 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
399 & ((blower(k,l,1),l=1,k),k=1,3)
401 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
402 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
403 & ((blower(k,l,j),l=1,k),k=1,3)
410 akl=akl+blower(k,m,j)*blower(l,m,j)
425 C Read torsional parameters in old format
427 read (itorp,*,end=113,err=113) ntortyp,nterm_old
428 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
429 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
434 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
440 write (iout,'(/a/)') 'Torsional constants:'
443 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
444 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
450 C Read torsional parameters
452 read (itorp,*,end=113,err=113) ntortyp
453 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
456 itortyp(i)=-itortyp(-1)
458 c write (iout,*) 'ntortyp',ntortyp
459 do i=0,ntortyp,ntortyp-1
460 do j=-ntortyp,ntortyp
461 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
463 nterm(-i,-j,iblock)=nterm(i,j,iblock)
464 nlor(-i,-j,iblock)=nlor(i,j,iblock)
467 do k=1,nterm(i,j,iblock)
468 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
470 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
471 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
472 v0ij=v0ij+si*v1(k,i,j,iblock)
475 do k=1,nlor(i,j,iblock)
476 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
477 & vlor2(k,i,j),vlor3(k,i,j)
478 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
485 write (iout,'(/a/)') 'Torsional constants:'
488 write (iout,*) 'ityp',i,' jtyp',j
489 write (iout,*) 'Fourier constants'
490 do k=1,nterm(i,j,iblock)
491 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
494 write (iout,*) 'Lorenz constants'
495 do k=1,nlor(i,j,iblock)
496 write (iout,'(3(1pe15.5))')
497 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
503 C 6/23/01 Read parameters for double torsionals
507 do j=-ntortyp+1,ntortyp-1
508 do k=-ntortyp+1,ntortyp-1
509 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
510 write (iout,*) "OK onelett",
513 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
514 & .or. t3.ne.toronelet(k)) then
515 write (iout,*) "Error in double torsional parameter file",
518 call MPI_Finalize(Ierror)
520 stop "Error in double torsional parameter file"
522 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
523 & ntermd_2(i,j,k,iblock)
524 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
525 & ntermd_1(i,j,k,iblock))
526 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
527 & ntermd_1(i,j,k,iblock))
528 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
529 & ntermd_1(i,j,k,iblock))
530 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
531 & ntermd_1(i,j,k,iblock))
532 C Matrix of D parameters for one dimesional foureir series
533 do l=1, ntermd_1(i,j,k,iblock)
534 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
535 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
536 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
537 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
539 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
540 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
541 & v2s(m,l,i,j,k,iblock),
542 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
543 C Matrix of D parameters for two dimesional fourier series
544 do l=1,ntermd_2(i,j,k,iblock)
546 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
547 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
548 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
549 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
558 write (iout,*) 'Constants for double torsionals'
560 do j=-ntortyp,ntortyp
561 do k=-ntortyp,ntortyp
562 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
563 & ' nsingle',ntermd_1(i,j,k,iblock),
564 & ' ndouble',ntermd_2(i,j,k,iblock)
566 write (iout,*) 'Single angles:'
567 do l=1,ntermd_1(i,j,k,iblock)
568 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
569 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
570 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
573 write (iout,*) 'Pairs of angles:'
574 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
575 do l=1,ntermd_2(i,j,k,iblock)
576 write (iout,'(i5,20f10.5)')
577 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
580 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
581 do l=1,ntermd_2(i,j,k,iblock)
582 write (iout,'(i5,20f10.5)')
583 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
584 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
593 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
594 C correlation energies.
596 read (isccor,*,end=119,err=119) nterm_sccor
601 read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
608 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
611 write (iout,*) 'ityp',i,' jtyp',j
613 write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
619 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
620 C interaction energy of the Gly, Ala, and Pro prototypes.
624 write (iout,*) "Coefficients of the cumulants"
626 read (ifourier,*) nloctyp
629 read (ifourier,*,end=115,err=115)
630 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
632 write (iout,*) 'Type',i
633 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
659 c Ctilde(1,1,i)=0.0d0
660 c Ctilde(1,2,i)=0.0d0
661 c Ctilde(2,1,i)=0.0d0
662 c Ctilde(2,2,i)=0.0d0
675 c Dtilde(1,1,i)=0.0d0
676 c Dtilde(1,2,i)=0.0d0
677 c Dtilde(2,1,i)=0.0d0
678 c Dtilde(2,2,i)=0.0d0
679 EE(1,1,i)= b(10)+b(11)
680 EE(2,2,i)=-b(10)+b(11)
681 EE(2,1,i)= b(12)-b(13)
682 EE(1,2,i)= b(12)+b(13)
687 c ee(2,1,i)=ee(1,2,i)
691 write (iout,*) 'Type',i
693 write(iout,*) B1(1,i),B1(2,i)
695 write(iout,*) B2(1,i),B2(2,i)
698 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
702 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
706 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
711 C Read electrostatic-interaction parameters
715 write (iout,'(/a)') 'Electrostatic interaction constants:'
716 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
717 & 'IT','JT','APP','BPP','AEL6','AEL3'
719 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
720 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
721 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
722 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
727 app (i,j)=epp(i,j)*rri*rri
728 bpp (i,j)=-2.0D0*epp(i,j)*rri
729 ael6(i,j)=elpp6(i,j)*4.2D0**6
730 ael3(i,j)=elpp3(i,j)*4.2D0**3
731 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
732 & ael6(i,j),ael3(i,j)
736 C Read side-chain interaction parameters.
738 read (isidep,*,end=117,err=117) ipot,expon
739 if (ipot.lt.1 .or. ipot.gt.5) then
740 write (iout,'(2a)') 'Error while reading SC interaction',
741 & 'potential file - unknown potential type.'
743 call MPI_Finalize(Ierror)
749 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
750 & ', exponents are ',expon,2*expon
751 goto (10,20,30,30,40) ipot
752 C----------------------- LJ potential ---------------------------------
753 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
754 & (sigma0(i),i=1,ntyp)
756 write (iout,'(/a/)') 'Parameters of the LJ potential:'
757 write (iout,'(a/)') 'The epsilon array:'
758 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
759 write (iout,'(/a)') 'One-body parameters:'
760 write (iout,'(a,4x,a)') 'residue','sigma'
761 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
764 C----------------------- LJK potential --------------------------------
765 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
766 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
768 write (iout,'(/a/)') 'Parameters of the LJK potential:'
769 write (iout,'(a/)') 'The epsilon array:'
770 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
771 write (iout,'(/a)') 'One-body parameters:'
772 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
773 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
777 C---------------------- GB or BP potential -----------------------------
778 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
779 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
781 C For the GB potential convert sigma'**2 into chi'
784 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
788 write (iout,'(/a/)') 'Parameters of the BP potential:'
789 write (iout,'(a/)') 'The epsilon array:'
790 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
791 write (iout,'(/a)') 'One-body parameters:'
792 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
794 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
795 & chip(i),alp(i),i=1,ntyp)
798 C--------------------- GBV potential -----------------------------------
799 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
800 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
801 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
803 write (iout,'(/a/)') 'Parameters of the GBV potential:'
804 write (iout,'(a/)') 'The epsilon array:'
805 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
806 write (iout,'(/a)') 'One-body parameters:'
807 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
808 & 's||/s_|_^2',' chip ',' alph '
809 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
810 & sigii(i),chip(i),alp(i),i=1,ntyp)
814 C-----------------------------------------------------------------------
815 C Calculate the "working" parameters of SC interactions.
823 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
824 sigma(j,i)=sigma(i,j)
825 rs0(i,j)=dwa16*sigma(i,j)
829 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
830 & 'Working parameters of the SC interactions:',
831 & ' a ',' b ',' augm ',' sigma ',' r0 ',
836 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
845 sigeps=dsign(1.0D0,epsij)
847 aa(i,j)=epsij*rrij*rrij
848 bb(i,j)=-sigeps*epsij*rrij
856 ratsig1=sigt2sq/sigt1sq
857 ratsig2=1.0D0/ratsig1
858 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
859 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
860 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
864 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
865 sigmaii(i,j)=rsum_max
866 sigmaii(j,i)=rsum_max
868 c sigmaii(i,j)=r0(i,j)
869 c sigmaii(j,i)=r0(i,j)
871 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
872 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
873 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
874 augm(i,j)=epsij*r_augm**(2*expon)
875 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
882 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
883 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
884 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
890 C Define the SC-p interaction constants (hard-coded; old style)
893 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
895 c aad(i,1)=0.3D0*4.0D0**12
896 C Following line for constants currently implemented
897 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
898 aad(i,1)=1.5D0*4.0D0**12
899 c aad(i,1)=0.17D0*5.6D0**12
901 C "Soft" SC-p repulsion
903 C Following line for constants currently implemented
904 c aad(i,1)=0.3D0*4.0D0**6
905 C "Hard" SC-p repulsion
906 bad(i,1)=3.0D0*4.0D0**6
907 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
916 C 8/9/01 Read the SC-p interaction constants from file
919 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
922 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
923 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
924 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
925 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
929 write (iout,*) "Parameters of SC-p interactions:"
931 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
932 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
937 C Define the constants of the disulfide bridge
941 c Old arbitrary potential - commented out.
946 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
947 c energy surface of diethyl disulfide.
948 c A. Liwo and U. Kozlowska, 11/24/03
965 write (iout,'(/a)') "Disulfide bridge parameters:"
966 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
967 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
968 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
969 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
973 111 write (iout,*) "Error reading bending energy parameters."
975 112 write (iout,*) "Error reading rotamer energy parameters."
977 113 write (iout,*) "Error reading torsional energy parameters."
979 114 write (iout,*) "Error reading double torsional energy parameters."
982 & "Error reading cumulant (multibody energy) parameters."
984 116 write (iout,*) "Error reading electrostatic energy parameters."
986 117 write (iout,*) "Error reading side chain interaction parameters."
988 118 write (iout,*) "Error reading SCp interaction parameters."
990 119 write (iout,*) "Error reading SCCOR parameters"
993 call MPI_Finalize(Ierror)
1000 subroutine getenv_loc(var, val)
1001 character(*) var, val
1004 character(2000) line
1007 open (196,file='env',status='old',readonly,shared)
1009 c write(*,*)'looking for ',var
1010 10 read(196,*,err=11,end=11)line
1011 iread=index(line,var)
1012 c write(*,*)iread,' ',var,' ',line
1013 if (iread.eq.0) go to 10
1014 c write(*,*)'---> ',line
1020 iread=iread+ilen(var)+1
1021 read (line(iread:),*,err=12,end=12) val
1022 c write(*,*)'OK: ',var,' = ',val
1028 #elif (defined CRAY)
1029 integer lennam,lenval,ierror
1031 c getenv using a POSIX call, useful on the T3D
1032 c Sept 1996, comment out error check on advice of H. Pritchard
1035 if(lennam.le.0) stop '--error calling getenv--'
1036 call pxfgetenv(var,lennam,val,lenval,ierror)
1037 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1039 call getenv(var,val)