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)
529 nterm_sccor(-i,j)=nterm_sccor(i,j)
530 nterm_sccor(-i,-j)=nterm_sccor(i,j)
531 nterm_sccor(i,-j)=nterm_sccor(i,j)
532 do k=1,nterm_sccor(i,j)
533 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
535 if (j.eq.iscprol) then
536 if (i.eq.isccortyp(10)) then
537 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
538 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
540 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
541 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
542 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
543 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
544 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
545 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
546 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
547 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
550 if (i.eq.isccortyp(10)) then
551 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
552 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
554 if (j.eq.isccortyp(10)) then
555 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
556 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
558 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
559 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
560 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
561 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
562 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
563 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
567 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
568 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
569 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
570 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
573 do k=1,nlor_sccor(i,j)
574 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
575 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
576 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
577 &(1+vlor3sccor(k,i,j)**2)
579 v0sccor(l,i,j)=v0ijsccor
580 v0sccor(l,-i,j)=v0ijsccor1
581 v0sccor(l,i,-j)=v0ijsccor2
582 v0sccor(l,-i,-j)=v0ijsccor3
588 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
589 c write (iout,*) 'ntortyp',ntortyp
591 cc maxinter is maximum interaction sites
595 read (isccor,*,end=113,err=113)
596 & nterm_sccor(i,j),nlor_sccor(i,j)
600 do k=1,nterm_sccor(i,j)
601 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
603 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
606 do k=1,nlor_sccor(i,j)
607 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
608 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
609 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
610 &(1+vlor3sccor(k,i,j)**2)
612 v0sccor(i,j)=v0ijsccor
620 write (iout,'(/a/)') 'Torsional constants:'
623 write (iout,*) 'ityp',i,' jtyp',j
624 write (iout,*) 'Fourier constants'
625 do k=1,nterm_sccor(i,j)
626 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
628 write (iout,*) 'Lorenz constants'
629 do k=1,nlor_sccor(i,j)
630 write (iout,'(3(1pe15.5))')
631 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
638 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
639 C interaction energy of the Gly, Ala, and Pro prototypes.
643 write (iout,*) "Coefficients of the cumulants"
645 read (ifourier,*) nloctyp
647 read (ifourier,*,end=115,err=115)
648 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
650 write (iout,*) 'Type',i
651 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
677 c Ctilde(1,1,i)=0.0d0
678 c Ctilde(1,2,i)=0.0d0
679 c Ctilde(2,1,i)=0.0d0
680 c Ctilde(2,2,i)=0.0d0
693 c Dtilde(1,1,i)=0.0d0
694 c Dtilde(1,2,i)=0.0d0
695 c Dtilde(2,1,i)=0.0d0
696 c Dtilde(2,2,i)=0.0d0
697 EE(1,1,i)= b(10)+b(11)
698 EE(2,2,i)=-b(10)+b(11)
699 EE(2,1,i)= b(12)-b(13)
700 EE(1,2,i)= b(12)+b(13)
705 c ee(2,1,i)=ee(1,2,i)
709 write (iout,*) 'Type',i
711 write(iout,*) B1(1,i),B1(2,i)
713 write(iout,*) B2(1,i),B2(2,i)
716 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
720 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
724 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
729 C Read electrostatic-interaction parameters
733 write (iout,'(/a)') 'Electrostatic interaction constants:'
734 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
735 & 'IT','JT','APP','BPP','AEL6','AEL3'
737 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
738 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
739 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
740 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
745 app (i,j)=epp(i,j)*rri*rri
746 bpp (i,j)=-2.0D0*epp(i,j)*rri
747 ael6(i,j)=elpp6(i,j)*4.2D0**6
748 ael3(i,j)=elpp3(i,j)*4.2D0**3
749 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
750 & ael6(i,j),ael3(i,j)
754 C Read side-chain interaction parameters.
756 read (isidep,*,end=117,err=117) ipot,expon
757 if (ipot.lt.1 .or. ipot.gt.5) then
758 write (iout,'(2a)') 'Error while reading SC interaction',
759 & 'potential file - unknown potential type.'
761 call MPI_Finalize(Ierror)
767 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
768 & ', exponents are ',expon,2*expon
769 goto (10,20,30,30,40) ipot
770 C----------------------- LJ potential ---------------------------------
771 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
772 & (sigma0(i),i=1,ntyp)
774 write (iout,'(/a/)') 'Parameters of the LJ potential:'
775 write (iout,'(a/)') 'The epsilon array:'
776 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
777 write (iout,'(/a)') 'One-body parameters:'
778 write (iout,'(a,4x,a)') 'residue','sigma'
779 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
782 C----------------------- LJK potential --------------------------------
783 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
784 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
786 write (iout,'(/a/)') 'Parameters of the LJK potential:'
787 write (iout,'(a/)') 'The epsilon array:'
788 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
789 write (iout,'(/a)') 'One-body parameters:'
790 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
791 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
795 C---------------------- GB or BP potential -----------------------------
796 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
797 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
799 C For the GB potential convert sigma'**2 into chi'
802 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
806 write (iout,'(/a/)') 'Parameters of the BP potential:'
807 write (iout,'(a/)') 'The epsilon array:'
808 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
809 write (iout,'(/a)') 'One-body parameters:'
810 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
812 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
813 & chip(i),alp(i),i=1,ntyp)
816 C--------------------- GBV potential -----------------------------------
817 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
818 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
819 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
821 write (iout,'(/a/)') 'Parameters of the GBV potential:'
822 write (iout,'(a/)') 'The epsilon array:'
823 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
824 write (iout,'(/a)') 'One-body parameters:'
825 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
826 & 's||/s_|_^2',' chip ',' alph '
827 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
828 & sigii(i),chip(i),alp(i),i=1,ntyp)
832 C-----------------------------------------------------------------------
833 C Calculate the "working" parameters of SC interactions.
841 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
842 sigma(j,i)=sigma(i,j)
843 rs0(i,j)=dwa16*sigma(i,j)
847 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
848 & 'Working parameters of the SC interactions:',
849 & ' a ',' b ',' augm ',' sigma ',' r0 ',
854 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
863 sigeps=dsign(1.0D0,epsij)
865 aa(i,j)=epsij*rrij*rrij
866 bb(i,j)=-sigeps*epsij*rrij
874 ratsig1=sigt2sq/sigt1sq
875 ratsig2=1.0D0/ratsig1
876 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
877 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
878 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
882 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
883 sigmaii(i,j)=rsum_max
884 sigmaii(j,i)=rsum_max
886 c sigmaii(i,j)=r0(i,j)
887 c sigmaii(j,i)=r0(i,j)
889 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
890 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
891 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
892 augm(i,j)=epsij*r_augm**(2*expon)
893 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
900 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
901 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
902 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
908 C Define the SC-p interaction constants (hard-coded; old style)
911 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
913 c aad(i,1)=0.3D0*4.0D0**12
914 C Following line for constants currently implemented
915 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
916 aad(i,1)=1.5D0*4.0D0**12
917 c aad(i,1)=0.17D0*5.6D0**12
919 C "Soft" SC-p repulsion
921 C Following line for constants currently implemented
922 c aad(i,1)=0.3D0*4.0D0**6
923 C "Hard" SC-p repulsion
924 bad(i,1)=3.0D0*4.0D0**6
925 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
934 C 8/9/01 Read the SC-p interaction constants from file
937 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
940 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
941 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
942 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
943 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
947 write (iout,*) "Parameters of SC-p interactions:"
949 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
950 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
955 C Define the constants of the disulfide bridge
959 c Old arbitrary potential - commented out.
964 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
965 c energy surface of diethyl disulfide.
966 c A. Liwo and U. Kozlowska, 11/24/03
983 write (iout,'(/a)') "Disulfide bridge parameters:"
984 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
985 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
986 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
987 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
991 111 write (iout,*) "Error reading bending energy parameters."
993 112 write (iout,*) "Error reading rotamer energy parameters."
995 113 write (iout,*) "Error reading torsional energy parameters."
997 114 write (iout,*) "Error reading double torsional energy parameters."
1000 & "Error reading cumulant (multibody energy) parameters."
1002 116 write (iout,*) "Error reading electrostatic energy parameters."
1004 117 write (iout,*) "Error reading side chain interaction parameters."
1006 118 write (iout,*) "Error reading SCp interaction parameters."
1008 119 write (iout,*) "Error reading SCCOR parameters"
1011 call MPI_Finalize(Ierror)
1018 subroutine getenv_loc(var, val)
1019 character(*) var, val
1022 character(2000) line
1025 open (196,file='env',status='old',readonly,shared)
1027 c write(*,*)'looking for ',var
1028 10 read(196,*,err=11,end=11)line
1029 iread=index(line,var)
1030 c write(*,*)iread,' ',var,' ',line
1031 if (iread.eq.0) go to 10
1032 c write(*,*)'---> ',line
1038 iread=iread+ilen(var)+1
1039 read (line(iread:),*,err=12,end=12) val
1040 c write(*,*)'OK: ',var,' = ',val
1046 #elif (defined CRAY)
1047 integer lennam,lenval,ierror
1049 c getenv using a POSIX call, useful on the T3D
1050 c Sept 1996, comment out error check on advice of H. Pritchard
1053 if(lennam.le.0) stop '--error calling getenv--'
1054 call pxfgetenv(var,lennam,val,lenval,ierror)
1055 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1057 call getenv(var,val)