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)
283 C Read the parameters of the probability distribution/energy expression
284 C of the side chains.
287 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
291 dsc_inv(i)=1.0D0/dsc(i)
302 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
303 & ((blower(k,l,1),l=1,k),k=1,3)
305 read (irotam,*,end=112,err=112) bsc(j,i)
306 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
307 & ((blower(k,l,j),l=1,k),k=1,3)
314 akl=akl+blower(k,m,j)*blower(l,m,j)
325 write (iout,'(/a)') 'Parameters of side-chain local geometry'
330 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
331 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
332 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
333 & 'log h',(bsc(j,i),j=1,nlobi)
334 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
335 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
337 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
338 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
341 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
342 write (iout,'(a,f10.4,4(16x,f10.4))')
343 & 'Center ',(bsc(j,i),j=1,nlobi)
344 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
353 C Read scrot parameters for potentials determined from all-atom AM1 calculations
354 C added by Urszula Kozlowska 07/11/2007
357 read (irotam,*,end=112,err=112)
359 read (irotam,*,end=112,err=112)
362 read(irotam,*,end=112,err=112) sc_parmin(j,i)
371 C Read torsional parameters in old format
373 read (itorp,*,end=113,err=113) ntortyp,nterm_old
374 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
375 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
380 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
386 write (iout,'(/a/)') 'Torsional constants:'
389 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
390 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
396 C Read torsional parameters
398 read (itorp,*,end=113,err=113) ntortyp
399 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
400 c write (iout,*) 'ntortyp',ntortyp
403 read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
407 read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
408 v0ij=v0ij+si*v1(k,i,j)
412 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
413 & vlor2(k,i,j),vlor3(k,i,j)
414 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
421 write (iout,'(/a/)') 'Torsional constants:'
424 write (iout,*) 'ityp',i,' jtyp',j
425 write (iout,*) 'Fourier constants'
427 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
429 write (iout,*) 'Lorenz constants'
431 write (iout,'(3(1pe15.5))')
432 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
438 C 6/23/01 Read parameters for double torsionals
443 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
444 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
445 & .or. t3.ne.onelett(k)) then
446 write (iout,*) "Error in double torsional parameter file",
449 call MPI_Finalize(Ierror)
451 stop "Error in double torsional parameter file"
453 read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
455 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
457 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
459 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
461 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
463 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
464 & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
465 & m=1,l-1),l=1,ntermd_2(i,j,k))
471 write (iout,*) 'Constants for double torsionals'
475 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
476 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
478 write (iout,*) 'Single angles:'
479 do l=1,ntermd_1(i,j,k)
480 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
481 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
482 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
485 write (iout,*) 'Pairs of angles:'
486 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
487 do l=1,ntermd_2(i,j,k)
488 write (iout,'(i5,20f10.5)')
489 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
492 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
493 do l=1,ntermd_2(i,j,k)
494 write (iout,'(i5,20f10.5)')
495 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
503 C Read of Side-chain backbone correlation parameters
504 C Modified 11 May 2012 by Adasko
507 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
509 read (isccor,*,end=119,err=119) nsccortyp
511 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
513 isccortyp(i)=-isccortyp(-i)
515 iscprol=isccortyp(20)
516 c write (iout,*) 'ntortyp',ntortyp
518 cc maxinter is maximum interaction sites
522 read (isccor,*,end=119,err=119)
523 &nterm_sccor(i,j),nlor_sccor(i,j)
526 nterm_sccor(-i,j)=nterm_sccor(i,j)
527 nterm_sccor(-i,-j)=nterm_sccor(i,j)
528 nterm_sccor(i,-j)=nterm_sccor(i,j)
529 do k=1,nterm_sccor(i,j)
530 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
532 if (j.eq.iscprol) then
533 if (i.eq.isccortyp(10)) then
534 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
535 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
537 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
538 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
539 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
540 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
541 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
542 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
543 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
544 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
547 if (i.eq.isccortyp(10)) then
548 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
549 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
551 if (j.eq.isccortyp(10)) then
552 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
553 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
555 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
556 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
557 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
558 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
559 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
560 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
564 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
567 do k=1,nlor_sccor(i,j)
568 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
569 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
570 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
571 &(1+vlor3sccor(k,i,j)**2)
573 v0sccor(i,j)=v0ijsccor
579 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
580 c write (iout,*) 'ntortyp',ntortyp
582 cc maxinter is maximum interaction sites
586 read (isccor,*,end=113,err=113)
587 & nterm_sccor(i,j),nlor_sccor(i,j)
591 do k=1,nterm_sccor(i,j)
592 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
594 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
597 do k=1,nlor_sccor(i,j)
598 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
599 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
600 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
601 &(1+vlor3sccor(k,i,j)**2)
603 v0sccor(i,j)=v0ijsccor
611 write (iout,'(/a/)') 'Torsional constants:'
614 write (iout,*) 'ityp',i,' jtyp',j
615 write (iout,*) 'Fourier constants'
616 do k=1,nterm_sccor(i,j)
617 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
619 write (iout,*) 'Lorenz constants'
620 do k=1,nlor_sccor(i,j)
621 write (iout,'(3(1pe15.5))')
622 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
629 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
630 C interaction energy of the Gly, Ala, and Pro prototypes.
634 write (iout,*) "Coefficients of the cumulants"
636 read (ifourier,*) nloctyp
638 read (ifourier,*,end=115,err=115)
639 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
641 write (iout,*) 'Type',i
642 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
668 c Ctilde(1,1,i)=0.0d0
669 c Ctilde(1,2,i)=0.0d0
670 c Ctilde(2,1,i)=0.0d0
671 c Ctilde(2,2,i)=0.0d0
684 c Dtilde(1,1,i)=0.0d0
685 c Dtilde(1,2,i)=0.0d0
686 c Dtilde(2,1,i)=0.0d0
687 c Dtilde(2,2,i)=0.0d0
688 EE(1,1,i)= b(10)+b(11)
689 EE(2,2,i)=-b(10)+b(11)
690 EE(2,1,i)= b(12)-b(13)
691 EE(1,2,i)= b(12)+b(13)
696 c ee(2,1,i)=ee(1,2,i)
700 write (iout,*) 'Type',i
702 write(iout,*) B1(1,i),B1(2,i)
704 write(iout,*) B2(1,i),B2(2,i)
707 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
711 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
715 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
720 C Read electrostatic-interaction parameters
724 write (iout,'(/a)') 'Electrostatic interaction constants:'
725 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
726 & 'IT','JT','APP','BPP','AEL6','AEL3'
728 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
729 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
730 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
731 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
736 app (i,j)=epp(i,j)*rri*rri
737 bpp (i,j)=-2.0D0*epp(i,j)*rri
738 ael6(i,j)=elpp6(i,j)*4.2D0**6
739 ael3(i,j)=elpp3(i,j)*4.2D0**3
740 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
741 & ael6(i,j),ael3(i,j)
745 C Read side-chain interaction parameters.
747 read (isidep,*,end=117,err=117) ipot,expon
748 if (ipot.lt.1 .or. ipot.gt.5) then
749 write (iout,'(2a)') 'Error while reading SC interaction',
750 & 'potential file - unknown potential type.'
752 call MPI_Finalize(Ierror)
758 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
759 & ', exponents are ',expon,2*expon
760 goto (10,20,30,30,40) ipot
761 C----------------------- LJ potential ---------------------------------
762 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
763 & (sigma0(i),i=1,ntyp)
765 write (iout,'(/a/)') 'Parameters of the LJ potential:'
766 write (iout,'(a/)') 'The epsilon array:'
767 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
768 write (iout,'(/a)') 'One-body parameters:'
769 write (iout,'(a,4x,a)') 'residue','sigma'
770 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
773 C----------------------- LJK potential --------------------------------
774 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
775 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
777 write (iout,'(/a/)') 'Parameters of the LJK potential:'
778 write (iout,'(a/)') 'The epsilon array:'
779 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
780 write (iout,'(/a)') 'One-body parameters:'
781 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
782 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
786 C---------------------- GB or BP potential -----------------------------
787 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
788 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
790 C For the GB potential convert sigma'**2 into chi'
793 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
797 write (iout,'(/a/)') 'Parameters of the BP potential:'
798 write (iout,'(a/)') 'The epsilon array:'
799 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
800 write (iout,'(/a)') 'One-body parameters:'
801 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
803 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
804 & chip(i),alp(i),i=1,ntyp)
807 C--------------------- GBV potential -----------------------------------
808 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
809 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
810 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
812 write (iout,'(/a/)') 'Parameters of the GBV potential:'
813 write (iout,'(a/)') 'The epsilon array:'
814 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
815 write (iout,'(/a)') 'One-body parameters:'
816 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
817 & 's||/s_|_^2',' chip ',' alph '
818 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
819 & sigii(i),chip(i),alp(i),i=1,ntyp)
823 C-----------------------------------------------------------------------
824 C Calculate the "working" parameters of SC interactions.
832 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
833 sigma(j,i)=sigma(i,j)
834 rs0(i,j)=dwa16*sigma(i,j)
838 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
839 & 'Working parameters of the SC interactions:',
840 & ' a ',' b ',' augm ',' sigma ',' r0 ',
845 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
854 sigeps=dsign(1.0D0,epsij)
856 aa(i,j)=epsij*rrij*rrij
857 bb(i,j)=-sigeps*epsij*rrij
865 ratsig1=sigt2sq/sigt1sq
866 ratsig2=1.0D0/ratsig1
867 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
868 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
869 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
873 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
874 sigmaii(i,j)=rsum_max
875 sigmaii(j,i)=rsum_max
877 c sigmaii(i,j)=r0(i,j)
878 c sigmaii(j,i)=r0(i,j)
880 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
881 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
882 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
883 augm(i,j)=epsij*r_augm**(2*expon)
884 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
891 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
892 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
893 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
899 C Define the SC-p interaction constants (hard-coded; old style)
902 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
904 c aad(i,1)=0.3D0*4.0D0**12
905 C Following line for constants currently implemented
906 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
907 aad(i,1)=1.5D0*4.0D0**12
908 c aad(i,1)=0.17D0*5.6D0**12
910 C "Soft" SC-p repulsion
912 C Following line for constants currently implemented
913 c aad(i,1)=0.3D0*4.0D0**6
914 C "Hard" SC-p repulsion
915 bad(i,1)=3.0D0*4.0D0**6
916 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
925 C 8/9/01 Read the SC-p interaction constants from file
928 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
931 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
932 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
933 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
934 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
938 write (iout,*) "Parameters of SC-p interactions:"
940 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
941 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
946 C Define the constants of the disulfide bridge
950 c Old arbitrary potential - commented out.
955 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
956 c energy surface of diethyl disulfide.
957 c A. Liwo and U. Kozlowska, 11/24/03
974 write (iout,'(/a)') "Disulfide bridge parameters:"
975 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
976 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
977 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
978 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
982 111 write (iout,*) "Error reading bending energy parameters."
984 112 write (iout,*) "Error reading rotamer energy parameters."
986 113 write (iout,*) "Error reading torsional energy parameters."
988 114 write (iout,*) "Error reading double torsional energy parameters."
991 & "Error reading cumulant (multibody energy) parameters."
993 116 write (iout,*) "Error reading electrostatic energy parameters."
995 117 write (iout,*) "Error reading side chain interaction parameters."
997 118 write (iout,*) "Error reading SCp interaction parameters."
999 119 write (iout,*) "Error reading SCCOR parameters"
1002 call MPI_Finalize(Ierror)
1009 subroutine getenv_loc(var, val)
1010 character(*) var, val
1013 character(2000) line
1016 open (196,file='env',status='old',readonly,shared)
1018 c write(*,*)'looking for ',var
1019 10 read(196,*,err=11,end=11)line
1020 iread=index(line,var)
1021 c write(*,*)iread,' ',var,' ',line
1022 if (iread.eq.0) go to 10
1023 c write(*,*)'---> ',line
1029 iread=iread+ilen(var)+1
1030 read (line(iread:),*,err=12,end=12) val
1031 c write(*,*)'OK: ',var,' = ',val
1037 #elif (defined CRAY)
1038 integer lennam,lenval,ierror
1040 c getenv using a POSIX call, useful on the T3D
1041 c Sept 1996, comment out error check on advice of H. Pritchard
1044 if(lennam.le.0) stop '--error calling getenv--'
1045 call pxfgetenv(var,lennam,val,lenval,ierror)
1046 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1048 call getenv(var,val)