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 implicit real*8 (a-h,o-z)
9 include 'COMMON.IOUNITS'
10 include 'COMMON.CHAIN'
11 include 'COMMON.INTERACT'
13 include 'COMMON.LOCAL'
14 include 'COMMON.TORSION'
15 include 'COMMON.FFIELD'
16 include 'COMMON.NAMES'
17 include 'COMMON.SBRIDGE'
18 include 'COMMON.SCCOR'
19 include 'COMMON.SCROT'
21 character*1 onelett(-2:2) /"p","a","G","A","P"/
23 dimension blower(3,3,maxlob)
24 double precision ip,mp
28 C Set LPRINT=.TRUE. for debugging
29 dwa16=2.0d0**(1.0d0/6.0d0)
32 C Assign virtual-bond length
37 read (ibond,*) vbldp0,akp
40 read (ibond,*) vbldsc0(1,i),aksc(1,i)
45 dsc_inv(i)=1.0D0/dsc(i)
49 read (ibond,*) ijunk,vbldp0,akp,rjunk
51 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
57 dsc_inv(i)=1.0D0/dsc(i)
62 write(iout,'(/a/)')"Force constants virtual bonds:"
63 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
65 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
67 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
68 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
70 write (iout,'(13x,3f10.5)')
71 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
77 C Read the parameters of the probability distribution/energy expression
78 C of the virtual-bond valence angles theta
81 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
82 & (bthet(j,i,1,1),j=1,2)
83 read (ithep,*) (polthet(j,i),j=0,3)
84 read (ithep,*) (gthet(j,i),j=1,3)
85 read (ithep,*) theta0(i),sig0(i),sigc0(i)
89 athet(1,i,1,-1)=athet(1,i,1,1)
90 athet(2,i,1,-1)=athet(2,i,1,1)
91 bthet(1,i,1,-1)=-bthet(1,i,1,1)
92 bthet(2,i,1,-1)=-bthet(2,i,1,1)
93 athet(1,i,-1,1)=-athet(1,i,1,1)
94 athet(2,i,-1,1)=-athet(2,i,1,1)
95 bthet(1,i,-1,1)=bthet(1,i,1,1)
96 bthet(2,i,-1,1)=bthet(2,i,1,1)
100 athet(1,i,-1,-1)=athet(1,-i,1,1)
101 athet(2,i,-1,-1)=-athet(2,-i,1,1)
102 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
103 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
104 athet(1,i,-1,1)=athet(1,-i,1,1)
105 athet(2,i,-1,1)=-athet(2,-i,1,1)
106 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
107 bthet(2,i,-1,1)=bthet(2,-i,1,1)
108 athet(1,i,1,-1)=-athet(1,-i,1,1)
109 athet(2,i,1,-1)=athet(2,-i,1,1)
110 bthet(1,i,1,-1)=bthet(1,-i,1,1)
111 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
116 polthet(j,i)=polthet(j,-i)
119 gthet(j,i)=gthet(j,-i)
125 c & 'Parameters of the virtual-bond valence angles:'
126 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
127 c & ' ATHETA0 ',' A1 ',' A2 ',
130 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
131 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
133 c write (iout,'(/a/9x,5a/79(1h-))')
134 c & 'Parameters of the expression for sigma(theta_c):',
135 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
136 c & ' ALPH3 ',' SIGMA0C '
138 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
139 c & (polthet(j,i),j=0,3),sigc0(i)
141 c write (iout,'(/a/9x,5a/79(1h-))')
142 c & 'Parameters of the second gaussian:',
143 c & ' THETA0 ',' SIGMA0 ',' G1 ',
146 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
147 c & sig0(i),(gthet(j,i),j=1,3)
150 & 'Parameters of the virtual-bond valence angles:'
151 write (iout,'(/a/9x,5a/79(1h-))')
152 & 'Coefficients of expansion',
153 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
154 & ' b1*10^1 ',' b2*10^1 '
156 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
157 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
158 & (10*bthet(j,i,1,1),j=1,2)
160 write (iout,'(/a/9x,5a/79(1h-))')
161 & 'Parameters of the expression for sigma(theta_c):',
162 & ' alpha0 ',' alph1 ',' alph2 ',
163 & ' alhp3 ',' sigma0c '
165 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
166 & (polthet(j,i),j=0,3),sigc0(i)
168 write (iout,'(/a/9x,5a/79(1h-))')
169 & 'Parameters of the second gaussian:',
170 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
173 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
174 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
179 C Read the parameters of Utheta determined from ab initio surfaces
180 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
182 read (ithep,*) nthetyp,ntheterm,ntheterm2,
183 & ntheterm3,nsingle,ndouble
184 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
185 read (ithep,*) (ithetyp(i),i=1,ntyp1)
187 ithetyp(i)=-ithetyp(-i)
190 do i=-maxthetyp,maxthetyp
191 do j=-maxthetyp,maxthetyp
192 do k=-maxthetyp,maxthetyp
193 aa0thet(i,j,k,iblock)=0.0d0
195 aathet(l,i,j,k,iblock)=0.0d0
199 bbthet(m,l,i,j,k,iblock)=0.0d0
200 ccthet(m,l,i,j,k,iblock)=0.0d0
201 ddthet(m,l,i,j,k,iblock)=0.0d0
202 eethet(m,l,i,j,k,iblock)=0.0d0
208 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
209 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
219 do j=-nthetyp,nthetyp
220 do k=-nthetyp,nthetyp
221 read (ithep,'(6a)') res1
222 read (ithep,*) aa0thet(i,j,k,iblock)
223 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
225 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
226 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
227 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
228 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
231 & (((ffthet(llll,lll,ll,i,j,k,iblock),
232 & ffthet(lll,llll,ll,i,j,k,iblock),
233 & ggthet(llll,lll,ll,i,j,k,iblock),
234 & ggthet(lll,llll,ll,i,j,k,iblock),
235 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
240 C For dummy ends assign glycine-type coefficients of theta-only terms; the
241 C coefficients of theta-and-gamma-dependent terms are zero.
246 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
247 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
249 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
250 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
253 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
255 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
258 C Substitution for D aminoacids from symmetry.
261 do j=-nthetyp,nthetyp
262 do k=-nthetyp,nthetyp
263 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
265 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
269 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
270 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
271 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
272 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
278 ffthet(llll,lll,ll,i,j,k,iblock)=
279 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
280 ffthet(lll,llll,ll,i,j,k,iblock)=
281 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
282 ggthet(llll,lll,ll,i,j,k,iblock)=
283 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
284 ggthet(lll,llll,ll,i,j,k,iblock)=
285 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
295 C Control printout of the coefficients of virtual-bond-angle potentials
298 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
302 write (iout,'(//4a)')
303 & 'Type ',onelett(i),onelett(j),onelett(k)
304 write (iout,'(//a,10x,a)') " l","a[l]"
305 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
306 write (iout,'(i2,1pe15.5)')
307 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
309 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
310 & "b",l,"c",l,"d",l,"e",l
312 write (iout,'(i2,4(1pe15.5))') m,
313 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
314 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
318 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
319 & "f+",l,"f-",l,"g+",l,"g-",l
322 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
323 & ffthet(n,m,l,i,j,k,iblock),
324 & ffthet(m,n,l,i,j,k,iblock),
325 & ggthet(n,m,l,i,j,k,iblock),
326 & ggthet(m,n,l,i,j,k,iblock)
339 C Read the parameters of the probability distribution/energy expression
340 C of the side chains.
343 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
347 dsc_inv(i)=1.0D0/dsc(i)
358 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
359 censc(1,1,-i)=censc(1,1,i)
360 censc(2,1,-i)=censc(2,1,i)
361 censc(3,1,-i)=-censc(3,1,i)
363 read (irotam,*) bsc(j,i)
364 read (irotam,*) (censc(k,j,i),k=1,3),
365 & ((blower(k,l,j),l=1,k),k=1,3)
366 censc(1,j,-i)=censc(1,j,i)
367 censc(2,j,-i)=censc(2,j,i)
368 censc(3,j,-i)=-censc(3,j,i)
369 C BSC is amplitude of Gaussian
376 akl=akl+blower(k,m,j)*blower(l,m,j)
380 if (((k.eq.3).and.(l.ne.3))
381 & .or.((l.eq.3).and.(k.ne.3))) then
382 gaussc(k,l,j,-i)=-akl
383 gaussc(l,k,j,-i)=-akl
395 write (iout,'(/a)') 'Parameters of side-chain local geometry'
399 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
400 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
401 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
402 c write (iout,'(a,f10.4,4(16x,f10.4))')
403 c & 'Center ',(bsc(j,i),j=1,nlobi)
404 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
405 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
406 & 'log h',(bsc(j,i),j=1,nlobi)
407 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
408 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
415 c blower(k,l,j)=gaussc(ind,j,i)
420 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
421 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
428 C Read scrot parameters for potentials determined from all-atom AM1 calculations
429 C added by Urszula Kozlowska 07/11/2007
437 read(irotam,*) sc_parmin(j,i)
445 C Read torsional parameters in old format
447 read (itorp,*) ntortyp,nterm_old
448 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
449 read (itorp,*) (itortyp(i),i=1,ntyp)
454 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
460 write (iout,'(/a/)') 'Torsional constants:'
463 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
464 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
470 C Read torsional parameters
472 read (itorp,*) ntortyp
473 read (itorp,*) (itortyp(i),i=1,ntyp)
474 write (iout,*) 'ntortyp',ntortyp
477 itortyp(i)=-itortyp(-i)
479 c write (iout,*) 'ntortyp',ntortyp
481 do j=-ntortyp+1,ntortyp-1
482 read (itorp,*) nterm(i,j,iblock),
484 nterm(-i,-j,iblock)=nterm(i,j,iblock)
485 nlor(-i,-j,iblock)=nlor(i,j,iblock)
488 do k=1,nterm(i,j,iblock)
489 read (itorp,*) kk,v1(k,i,j,iblock),
491 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
492 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
493 v0ij=v0ij+si*v1(k,i,j,iblock)
496 do k=1,nlor(i,j,iblock)
497 read (itorp,*) kk,vlor1(k,i,j),
498 & vlor2(k,i,j),vlor3(k,i,j)
499 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
502 v0(-i,-j,iblock)=v0ij
508 write (iout,'(/a/)') 'Torsional constants:'
511 write (iout,*) 'ityp',i,' jtyp',j
512 write (iout,*) 'Fourier constants'
513 do k=1,nterm(i,j,iblock)
514 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
517 write (iout,*) 'Lorenz constants'
518 do k=1,nlor(i,j,iblock)
519 write (iout,'(3(1pe15.5))')
520 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
526 C 6/23/01 Read parameters for double torsionals
530 do j=-ntortyp+1,ntortyp-1
531 do k=-ntortyp+1,ntortyp-1
532 read (itordp,'(3a1)') t1,t2,t3
533 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
534 & .or. t3.ne.onelett(k)) then
535 write (iout,*) "Error in double torsional parameter file",
537 stop "Error in double torsional parameter file"
539 read (itordp,*) ntermd_1(i,j,k,iblock),
540 & ntermd_2(i,j,k,iblock)
541 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
542 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
543 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
544 & ntermd_1(i,j,k,iblock))
545 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
546 & ntermd_1(i,j,k,iblock))
547 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
548 & ntermd_1(i,j,k,iblock))
549 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
550 & ntermd_1(i,j,k,iblock))
551 C Martix of D parameters for one dimesional foureir series
552 do l=1,ntermd_1(i,j,k,iblock)
553 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
554 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
555 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
556 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
557 c write(iout,*) "whcodze" ,
558 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
560 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
561 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
562 & v2s(m,l,i,j,k,iblock),
563 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
564 C Martix of D parameters for two dimesional fourier series
565 do l=1,ntermd_2(i,j,k,iblock)
567 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
568 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
569 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
570 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
579 write (iout,*) 'Constants for double torsionals'
582 do j=-ntortyp+1,ntortyp-1
583 do k=-ntortyp+1,ntortyp-1
584 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
585 & ' nsingle',ntermd_1(i,j,k,iblock),
586 & ' ndouble',ntermd_2(i,j,k,iblock)
588 write (iout,*) 'Single angles:'
589 do l=1,ntermd_1(i,j,k,iblock)
590 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
591 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
592 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
593 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
596 write (iout,*) 'Pairs of angles:'
597 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
598 do l=1,ntermd_2(i,j,k,iblock)
599 write (iout,'(i5,20f10.5)')
600 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
603 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
604 do l=1,ntermd_2(i,j,k,iblock)
605 write (iout,'(i5,20f10.5)')
606 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
607 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
616 C Read of Side-chain backbone correlation parameters
617 C Modified 11 May 2012 by Adasko
620 read (isccor,*) nsccortyp
621 read (isccor,*) (isccortyp(i),i=1,ntyp)
622 c write (iout,*) 'ntortyp',ntortyp
624 isccortyp(i)=-isccortyp(-i)
626 iscprol=isccortyp(20)
628 cc maxinter is maximum interaction sites
633 &nterm_sccor(i,j),nlor_sccor(i,j)
639 nterm_sccor(-i,j)=nterm_sccor(i,j)
640 nterm_sccor(-i,-j)=nterm_sccor(i,j)
641 nterm_sccor(i,-j)=nterm_sccor(i,j)
642 do k=1,nterm_sccor(i,j)
643 read (isccor,*) kk,v1sccor(k,l,i,j)
645 if (j.eq.iscprol) then
646 if (i.eq.isccortyp(10)) then
647 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
648 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
650 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
651 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
652 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
653 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
654 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
655 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
656 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
657 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
660 if (i.eq.isccortyp(10)) then
661 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
662 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
664 if (j.eq.isccortyp(10)) then
665 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
666 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
668 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
669 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
670 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
671 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
672 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
673 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
677 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
678 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
679 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
680 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
683 do k=1,nlor_sccor(i,j)
684 read (isccor,*) kk,vlor1sccor(k,i,j),
685 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
686 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
687 &(1+vlor3sccor(k,i,j)**2)
689 v0sccor(l,i,j)=v0ijsccor
690 v0sccor(l,-i,j)=v0ijsccor1
691 v0sccor(l,i,-j)=v0ijsccor2
692 v0sccor(l,-i,-j)=v0ijsccor3
699 write (iout,'(/a/)') 'Torsional constants:'
702 write (iout,*) 'ityp',i,' jtyp',j
703 write (iout,*) 'Fourier constants'
704 do k=1,nterm_sccor(i,j)
705 write (iout,'(2(1pe15.5))')
706 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
708 write (iout,*) 'Lorenz constants'
709 do k=1,nlor_sccor(i,j)
710 write (iout,'(3(1pe15.5))')
711 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
718 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
719 C interaction energy of the Gly, Ala, and Pro prototypes.
721 read (ifourier,*) nloctyp
724 read (ifourier,*) (b(ii,i),ii=1,13)
726 write (iout,*) 'Type',i
727 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
733 B1tilde(1,i) = b(3,i)
734 B1tilde(2,i) =-b(5,i)
735 B1tilde(1,-i) =-b(3,i)
736 B1tilde(2,-i) =b(5,i)
751 Ctilde(2,1,i)=-b(9,i)
753 Ctilde(1,1,-i)=b(7,i)
754 Ctilde(1,2,-i)=-b(9,i)
755 Ctilde(2,1,-i)=b(9,i)
756 Ctilde(2,2,-i)=b(7,i)
767 Dtilde(2,1,i)=-b(8,i)
769 Dtilde(1,1,-i)=b(6,i)
770 Dtilde(1,2,-i)=-b(8,i)
771 Dtilde(2,1,-i)=b(8,i)
772 Dtilde(2,2,-i)=b(6,i)
773 EE(1,1,i)= b(10,i)+b(11,i)
774 EE(2,2,i)=-b(10,i)+b(11,i)
775 EE(2,1,i)= b(12,i)-b(13,i)
776 EE(1,2,i)= b(12,i)+b(13,i)
777 EE(1,1,-i)= b(10,i)+b(11,i)
778 EE(2,2,-i)=-b(10,i)+b(11,i)
779 EE(2,1,-i)=-b(12,i)+b(13,i)
780 EE(1,2,-i)=-b(12,i)-b(13,i)
784 write (iout,*) 'Type',i
786 c write (iout,'(f10.5)') B1(:,i)
787 write(iout,*) B1(1,i),B1(2,i)
789 c write (iout,'(f10.5)') B2(:,i)
790 write(iout,*) B2(1,i),B2(2,i)
793 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
797 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
801 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
806 C Read electrostatic-interaction parameters
809 write (iout,'(/a)') 'Electrostatic interaction constants:'
810 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
811 & 'IT','JT','APP','BPP','AEL6','AEL3'
813 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
814 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
815 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
816 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
821 app (i,j)=epp(i,j)*rri*rri
822 bpp (i,j)=-2.0D0*epp(i,j)*rri
823 ael6(i,j)=elpp6(i,j)*4.2D0**6
824 ael3(i,j)=elpp3(i,j)*4.2D0**3
825 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
826 & ael6(i,j),ael3(i,j)
830 C Read side-chain interaction parameters.
832 read (isidep,*) ipot,expon
833 if (ipot.lt.1 .or. ipot.gt.5) then
834 write (iout,'(2a)') 'Error while reading SC interaction',
835 & 'potential file - unknown potential type.'
839 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
840 & ', exponents are ',expon,2*expon
841 goto (10,20,30,30,40) ipot
842 C----------------------- LJ potential ---------------------------------
843 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
845 write (iout,'(/a/)') 'Parameters of the LJ potential:'
846 write (iout,'(a/)') 'The epsilon array:'
847 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
848 write (iout,'(/a)') 'One-body parameters:'
849 write (iout,'(a,4x,a)') 'residue','sigma'
850 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
853 C----------------------- LJK potential --------------------------------
854 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
855 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
857 write (iout,'(/a/)') 'Parameters of the LJK potential:'
858 write (iout,'(a/)') 'The epsilon array:'
859 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
860 write (iout,'(/a)') 'One-body parameters:'
861 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
862 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
866 C---------------------- GB or BP potential -----------------------------
867 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
868 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
870 C For the GB potential convert sigma'**2 into chi'
873 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
877 write (iout,'(/a/)') 'Parameters of the BP potential:'
878 write (iout,'(a/)') 'The epsilon array:'
879 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
880 write (iout,'(/a)') 'One-body parameters:'
881 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
883 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
884 & chip(i),alp(i),i=1,ntyp)
887 C--------------------- GBV potential -----------------------------------
888 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
889 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
890 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
892 write (iout,'(/a/)') 'Parameters of the GBV potential:'
893 write (iout,'(a/)') 'The epsilon array:'
894 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
895 write (iout,'(/a)') 'One-body parameters:'
896 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
897 & 's||/s_|_^2',' chip ',' alph '
898 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
899 & sigii(i),chip(i),alp(i),i=1,ntyp)
903 C-----------------------------------------------------------------------
904 C Calculate the "working" parameters of SC interactions.
912 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
913 sigma(j,i)=sigma(i,j)
914 rs0(i,j)=dwa16*sigma(i,j)
918 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
919 & 'Working parameters of the SC interactions:',
920 & ' a ',' b ',' augm ',' sigma ',' r0 ',
925 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
934 sigeps=dsign(1.0D0,epsij)
936 aa(i,j)=epsij*rrij*rrij
937 bb(i,j)=-sigeps*epsij*rrij
945 ratsig1=sigt2sq/sigt1sq
946 ratsig2=1.0D0/ratsig1
947 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
948 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
949 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
953 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
954 sigmaii(i,j)=rsum_max
955 sigmaii(j,i)=rsum_max
957 c sigmaii(i,j)=r0(i,j)
958 c sigmaii(j,i)=r0(i,j)
960 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
961 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
962 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
963 augm(i,j)=epsij*r_augm**(2*expon)
964 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
971 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
972 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
973 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
978 C Define the SC-p interaction constants
988 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
990 c aad(i,1)=0.3D0*4.0D0**12
991 C Following line for constants currently implemented
992 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
993 aad(i,1)=1.5D0*4.0D0**12
994 c aad(i,1)=0.17D0*5.6D0**12
996 C "Soft" SC-p repulsion
998 C Following line for constants currently implemented
999 c aad(i,1)=0.3D0*4.0D0**6
1000 C "Hard" SC-p repulsion
1001 bad(i,1)=3.0D0*4.0D0**6
1002 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1011 C 8/9/01 Read the SC-p interaction constants from file
1014 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1017 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1018 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1019 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1020 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1024 write (iout,*) "Parameters of SC-p interactions:"
1026 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1027 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1032 C Define the constants of the disulfide bridge
1036 c Old arbitrary potential - commented out.
1041 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1042 c energy surface of diethyl disulfide.
1043 c A. Liwo and U. Kozlowska, 11/24/03