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)
191 aathet(l,i,j,k)=0.0d0
195 bbthet(m,l,i,j,k)=0.0d0
196 ccthet(m,l,i,j,k)=0.0d0
197 ddthet(m,l,i,j,k)=0.0d0
198 eethet(m,l,i,j,k)=0.0d0
204 ffthet(mm,m,l,i,j,k)=0.0d0
205 ggthet(mm,m,l,i,j,k)=0.0d0
215 read (ithep,'(3a)') res1,res2,res3
216 read (ithep,*) aa0thet(i,j,k)
217 read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
219 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
220 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
221 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
222 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
224 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
225 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
226 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
231 C For dummy ends assign glycine-type coefficients of theta-only terms; the
232 C coefficients of theta-and-gamma-dependent terms are zero.
237 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
238 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
240 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
241 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
244 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
246 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
249 C Control printout of the coefficients of virtual-bond-angle potentials
252 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
256 write (iout,'(//4a)')
257 & 'Type ',onelett(i),onelett(j),onelett(k)
258 write (iout,'(//a,10x,a)') " l","a[l]"
259 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
260 write (iout,'(i2,1pe15.5)')
261 & (l,aathet(l,i,j,k),l=1,ntheterm)
263 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
264 & "b",l,"c",l,"d",l,"e",l
266 write (iout,'(i2,4(1pe15.5))') m,
267 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
268 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
272 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
273 & "f+",l,"f-",l,"g+",l,"g-",l
276 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
277 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
278 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
291 C Read the parameters of the probability distribution/energy expression
292 C of the side chains.
295 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
299 dsc_inv(i)=1.0D0/dsc(i)
310 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
311 censc(1,1,-i)=censc(1,1,i)
312 censc(2,1,-i)=censc(2,1,i)
313 censc(3,1,-i)=-censc(3,1,i)
315 read (irotam,*) bsc(j,i)
316 read (irotam,*) (censc(k,j,i),k=1,3),
317 & ((blower(k,l,j),l=1,k),k=1,3)
318 censc(1,j,-i)=censc(1,j,i)
319 censc(2,j,-i)=censc(2,j,i)
320 censc(3,j,-i)=-censc(3,j,i)
321 C BSC is amplitude of Gaussian
328 akl=akl+blower(k,m,j)*blower(l,m,j)
332 if (((k.eq.3).and.(l.ne.3))
333 & .or.((l.eq.3).and.(k.ne.3))) then
334 gaussc(k,l,j,-i)=-akl
335 gaussc(l,k,j,-i)=-akl
347 write (iout,'(/a)') 'Parameters of side-chain local geometry'
351 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
352 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
353 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
354 c write (iout,'(a,f10.4,4(16x,f10.4))')
355 c & 'Center ',(bsc(j,i),j=1,nlobi)
356 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
357 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
358 & 'log h',(bsc(j,i),j=1,nlobi)
359 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
360 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
367 c blower(k,l,j)=gaussc(ind,j,i)
372 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
373 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
380 C Read scrot parameters for potentials determined from all-atom AM1 calculations
381 C added by Urszula Kozlowska 07/11/2007
389 read(irotam,*) sc_parmin(j,i)
397 C Read torsional parameters in old format
399 read (itorp,*) ntortyp,nterm_old
400 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
401 read (itorp,*) (itortyp(i),i=1,ntyp)
406 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
412 write (iout,'(/a/)') 'Torsional constants:'
415 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
416 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
422 C Read torsional parameters
424 read (itorp,*) ntortyp
425 read (itorp,*) (itortyp(i),i=1,ntyp)
426 write (iout,*) 'ntortyp',ntortyp
429 itortyp(i)=-itortyp(-i)
431 c write (iout,*) 'ntortyp',ntortyp
433 do j=-ntortyp+1,ntortyp-1
434 read (itorp,*) nterm(i,j,iblock),
436 nterm(-i,-j,iblock)=nterm(i,j,iblock)
437 nlor(-i,-j,iblock)=nlor(i,j,iblock)
440 do k=1,nterm(i,j,iblock)
441 read (itorp,*) kk,v1(k,i,j,iblock),
443 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
444 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
445 v0ij=v0ij+si*v1(k,i,j,iblock)
448 do k=1,nlor(i,j,iblock)
449 read (itorp,*) kk,vlor1(k,i,j),
450 & vlor2(k,i,j),vlor3(k,i,j)
451 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
454 v0(-i,-j,iblock)=v0ij
460 write (iout,'(/a/)') 'Torsional constants:'
463 write (iout,*) 'ityp',i,' jtyp',j
464 write (iout,*) 'Fourier constants'
465 do k=1,nterm(i,j,iblock)
466 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
469 write (iout,*) 'Lorenz constants'
470 do k=1,nlor(i,j,iblock)
471 write (iout,'(3(1pe15.5))')
472 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
478 C 6/23/01 Read parameters for double torsionals
482 do j=-ntortyp+1,ntortyp-1
483 do k=-ntortyp+1,ntortyp-1
484 read (itordp,'(3a1)') t1,t2,t3
485 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
486 & .or. t3.ne.onelett(k)) then
487 write (iout,*) "Error in double torsional parameter file",
489 stop "Error in double torsional parameter file"
491 read (itordp,*) ntermd_1(i,j,k,iblock),
492 & ntermd_2(i,j,k,iblock)
493 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
494 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
495 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
496 & ntermd_1(i,j,k,iblock))
497 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
498 & ntermd_1(i,j,k,iblock))
499 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
500 & ntermd_1(i,j,k,iblock))
501 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
502 & ntermd_1(i,j,k,iblock))
503 C Martix of D parameters for one dimesional foureir series
504 do l=1,ntermd_1(i,j,k,iblock)
505 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
506 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
507 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
508 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
509 c write(iout,*) "whcodze" ,
510 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
512 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
513 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
514 & v2s(m,l,i,j,k,iblock),
515 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
516 C Martix of D parameters for two dimesional fourier series
517 do l=1,ntermd_2(i,j,k,iblock)
519 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
520 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
521 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
522 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
531 write (iout,*) 'Constants for double torsionals'
534 do j=-ntortyp+1,ntortyp-1
535 do k=-ntortyp+1,ntortyp-1
536 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
537 & ' nsingle',ntermd_1(i,j,k,iblock),
538 & ' ndouble',ntermd_2(i,j,k,iblock)
540 write (iout,*) 'Single angles:'
541 do l=1,ntermd_1(i,j,k,iblock)
542 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
543 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
544 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
545 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
548 write (iout,*) 'Pairs of angles:'
549 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
550 do l=1,ntermd_2(i,j,k,iblock)
551 write (iout,'(i5,20f10.5)')
552 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
555 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
556 do l=1,ntermd_2(i,j,k,iblock)
557 write (iout,'(i5,20f10.5)')
558 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
559 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
568 C Read of Side-chain backbone correlation parameters
569 C Modified 11 May 2012 by Adasko
572 read (isccor,*) nsccortyp
573 read (isccor,*) (isccortyp(i),i=1,ntyp)
574 c write (iout,*) 'ntortyp',ntortyp
576 isccortyp(i)=-isccortyp(-i)
578 iscprol=isccortyp(20)
580 cc maxinter is maximum interaction sites
585 &nterm_sccor(i,j),nlor_sccor(i,j)
591 nterm_sccor(-i,j)=nterm_sccor(i,j)
592 nterm_sccor(-i,-j)=nterm_sccor(i,j)
593 nterm_sccor(i,-j)=nterm_sccor(i,j)
594 do k=1,nterm_sccor(i,j)
595 read (isccor,*) kk,v1sccor(k,l,i,j)
597 if (j.eq.iscprol) then
598 if (i.eq.isccortyp(10)) then
599 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
600 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
602 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
603 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
604 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
605 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
606 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
607 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
608 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
609 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
612 if (i.eq.isccortyp(10)) then
613 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
614 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
616 if (j.eq.isccortyp(10)) then
617 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
618 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
620 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
621 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
622 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
623 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
624 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
625 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
629 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
630 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
631 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
632 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
635 do k=1,nlor_sccor(i,j)
636 read (isccor,*) kk,vlor1sccor(k,i,j),
637 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
638 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
639 &(1+vlor3sccor(k,i,j)**2)
641 v0sccor(l,i,j)=v0ijsccor
642 v0sccor(l,-i,j)=v0ijsccor1
643 v0sccor(l,i,-j)=v0ijsccor2
644 v0sccor(l,-i,-j)=v0ijsccor3
651 write (iout,'(/a/)') 'Torsional constants:'
654 write (iout,*) 'ityp',i,' jtyp',j
655 write (iout,*) 'Fourier constants'
656 do k=1,nterm_sccor(i,j)
657 write (iout,'(2(1pe15.5))')
658 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
660 write (iout,*) 'Lorenz constants'
661 do k=1,nlor_sccor(i,j)
662 write (iout,'(3(1pe15.5))')
663 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
670 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
671 C interaction energy of the Gly, Ala, and Pro prototypes.
673 read (ifourier,*) nloctyp
676 read (ifourier,*) (b(ii,i),ii=1,13)
678 write (iout,*) 'Type',i
679 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
685 B1tilde(1,i) = b(3,i)
686 B1tilde(2,i) =-b(5,i)
687 B1tilde(1,-i) =-b(3,i)
688 B1tilde(2,-i) =b(5,i)
703 Ctilde(2,1,i)=-b(9,i)
705 Ctilde(1,1,-i)=b(7,i)
706 Ctilde(1,2,-i)=-b(9,i)
707 Ctilde(2,1,-i)=b(9,i)
708 Ctilde(2,2,-i)=b(7,i)
719 Dtilde(2,1,i)=-b(8,i)
721 Dtilde(1,1,-i)=b(6,i)
722 Dtilde(1,2,-i)=-b(8,i)
723 Dtilde(2,1,-i)=b(8,i)
724 Dtilde(2,2,-i)=b(6,i)
725 EE(1,1,i)= b(10,i)+b(11,i)
726 EE(2,2,i)=-b(10,i)+b(11,i)
727 EE(2,1,i)= b(12,i)-b(13,i)
728 EE(1,2,i)= b(12,i)+b(13,i)
729 EE(1,1,-i)= b(10,i)+b(11,i)
730 EE(2,2,-i)=-b(10,i)+b(11,i)
731 EE(2,1,-i)=-b(12,i)+b(13,i)
732 EE(1,2,-i)=-b(12,i)-b(13,i)
736 write (iout,*) 'Type',i
738 c write (iout,'(f10.5)') B1(:,i)
739 write(iout,*) B1(1,i),B1(2,i)
741 c write (iout,'(f10.5)') B2(:,i)
742 write(iout,*) B2(1,i),B2(2,i)
745 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
749 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
753 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
758 C Read electrostatic-interaction parameters
761 write (iout,'(/a)') 'Electrostatic interaction constants:'
762 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
763 & 'IT','JT','APP','BPP','AEL6','AEL3'
765 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
766 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
767 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
768 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
773 app (i,j)=epp(i,j)*rri*rri
774 bpp (i,j)=-2.0D0*epp(i,j)*rri
775 ael6(i,j)=elpp6(i,j)*4.2D0**6
776 ael3(i,j)=elpp3(i,j)*4.2D0**3
777 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
778 & ael6(i,j),ael3(i,j)
782 C Read side-chain interaction parameters.
784 read (isidep,*) ipot,expon
785 if (ipot.lt.1 .or. ipot.gt.5) then
786 write (iout,'(2a)') 'Error while reading SC interaction',
787 & 'potential file - unknown potential type.'
791 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
792 & ', exponents are ',expon,2*expon
793 goto (10,20,30,30,40) ipot
794 C----------------------- LJ potential ---------------------------------
795 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
797 write (iout,'(/a/)') 'Parameters of the LJ 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,a)') 'residue','sigma'
802 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
805 C----------------------- LJK potential --------------------------------
806 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
807 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
809 write (iout,'(/a/)') 'Parameters of the LJK potential:'
810 write (iout,'(a/)') 'The epsilon array:'
811 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
812 write (iout,'(/a)') 'One-body parameters:'
813 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
814 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
818 C---------------------- GB or BP potential -----------------------------
819 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
820 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
822 C For the GB potential convert sigma'**2 into chi'
825 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
829 write (iout,'(/a/)') 'Parameters of the BP potential:'
830 write (iout,'(a/)') 'The epsilon array:'
831 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
832 write (iout,'(/a)') 'One-body parameters:'
833 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
835 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
836 & chip(i),alp(i),i=1,ntyp)
839 C--------------------- GBV potential -----------------------------------
840 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
841 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
842 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
844 write (iout,'(/a/)') 'Parameters of the GBV potential:'
845 write (iout,'(a/)') 'The epsilon array:'
846 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
847 write (iout,'(/a)') 'One-body parameters:'
848 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
849 & 's||/s_|_^2',' chip ',' alph '
850 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
851 & sigii(i),chip(i),alp(i),i=1,ntyp)
855 C-----------------------------------------------------------------------
856 C Calculate the "working" parameters of SC interactions.
864 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
865 sigma(j,i)=sigma(i,j)
866 rs0(i,j)=dwa16*sigma(i,j)
870 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
871 & 'Working parameters of the SC interactions:',
872 & ' a ',' b ',' augm ',' sigma ',' r0 ',
877 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
886 sigeps=dsign(1.0D0,epsij)
888 aa(i,j)=epsij*rrij*rrij
889 bb(i,j)=-sigeps*epsij*rrij
897 ratsig1=sigt2sq/sigt1sq
898 ratsig2=1.0D0/ratsig1
899 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
900 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
901 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
905 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
906 sigmaii(i,j)=rsum_max
907 sigmaii(j,i)=rsum_max
909 c sigmaii(i,j)=r0(i,j)
910 c sigmaii(j,i)=r0(i,j)
912 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
913 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
914 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
915 augm(i,j)=epsij*r_augm**(2*expon)
916 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
923 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
924 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
925 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
930 C Define the SC-p interaction constants
940 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
942 c aad(i,1)=0.3D0*4.0D0**12
943 C Following line for constants currently implemented
944 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
945 aad(i,1)=1.5D0*4.0D0**12
946 c aad(i,1)=0.17D0*5.6D0**12
948 C "Soft" SC-p repulsion
950 C Following line for constants currently implemented
951 c aad(i,1)=0.3D0*4.0D0**6
952 C "Hard" SC-p repulsion
953 bad(i,1)=3.0D0*4.0D0**6
954 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
963 C 8/9/01 Read the SC-p interaction constants from file
966 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
969 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
970 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
971 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
972 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
976 write (iout,*) "Parameters of SC-p interactions:"
978 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
979 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
984 C Define the constants of the disulfide bridge
988 c Old arbitrary potential - commented out.
993 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
994 c energy surface of diethyl disulfide.
995 c A. Liwo and U. Kozlowska, 11/24/03
1005 write (iout,'(/a)') "Disulfide bridge parameters:"
1006 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1007 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1008 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1009 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,