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))
569 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
570 C correlation energies.
572 read (isccor,*) nterm_sccor
578 & kk,v1sccor(k,i,j),v2sccor(k,i,j)
584 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
587 write (iout,*) 'ityp',i,' jtyp',j
589 write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
595 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
596 C interaction energy of the Gly, Ala, and Pro prototypes.
598 read (ifourier,*) nloctyp
601 read (ifourier,*) (b(ii,i),ii=1,13)
603 write (iout,*) 'Type',i
604 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
610 B1tilde(1,i) = b(3,i)
611 B1tilde(2,i) =-b(5,i)
612 B1tilde(1,-i) =-b(3,i)
613 B1tilde(2,-i) =b(5,i)
628 Ctilde(2,1,i)=-b(9,i)
630 Ctilde(1,1,-i)=b(7,i)
631 Ctilde(1,2,-i)=-b(9,i)
632 Ctilde(2,1,-i)=b(9,i)
633 Ctilde(2,2,-i)=b(7,i)
644 Dtilde(2,1,i)=-b(8,i)
646 Dtilde(1,1,-i)=b(6,i)
647 Dtilde(1,2,-i)=-b(8,i)
648 Dtilde(2,1,-i)=b(8,i)
649 Dtilde(2,2,-i)=b(6,i)
650 EE(1,1,i)= b(10,i)+b(11,i)
651 EE(2,2,i)=-b(10,i)+b(11,i)
652 EE(2,1,i)= b(12,i)-b(13,i)
653 EE(1,2,i)= b(12,i)+b(13,i)
654 EE(1,1,-i)= b(10,i)+b(11,i)
655 EE(2,2,-i)=-b(10,i)+b(11,i)
656 EE(2,1,-i)=-b(12,i)+b(13,i)
657 EE(1,2,-i)=-b(12,i)-b(13,i)
661 write (iout,*) 'Type',i
663 c write (iout,'(f10.5)') B1(:,i)
664 write(iout,*) B1(1,i),B1(2,i)
666 c write (iout,'(f10.5)') B2(:,i)
667 write(iout,*) B2(1,i),B2(2,i)
670 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
674 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
678 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
683 C Read electrostatic-interaction parameters
686 write (iout,'(/a)') 'Electrostatic interaction constants:'
687 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
688 & 'IT','JT','APP','BPP','AEL6','AEL3'
690 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
691 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
692 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
693 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
698 app (i,j)=epp(i,j)*rri*rri
699 bpp (i,j)=-2.0D0*epp(i,j)*rri
700 ael6(i,j)=elpp6(i,j)*4.2D0**6
701 ael3(i,j)=elpp3(i,j)*4.2D0**3
702 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
703 & ael6(i,j),ael3(i,j)
707 C Read side-chain interaction parameters.
709 read (isidep,*) ipot,expon
710 if (ipot.lt.1 .or. ipot.gt.5) then
711 write (iout,'(2a)') 'Error while reading SC interaction',
712 & 'potential file - unknown potential type.'
716 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
717 & ', exponents are ',expon,2*expon
718 goto (10,20,30,30,40) ipot
719 C----------------------- LJ potential ---------------------------------
720 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
722 write (iout,'(/a/)') 'Parameters of the LJ potential:'
723 write (iout,'(a/)') 'The epsilon array:'
724 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
725 write (iout,'(/a)') 'One-body parameters:'
726 write (iout,'(a,4x,a)') 'residue','sigma'
727 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
730 C----------------------- LJK potential --------------------------------
731 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
732 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
734 write (iout,'(/a/)') 'Parameters of the LJK potential:'
735 write (iout,'(a/)') 'The epsilon array:'
736 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
737 write (iout,'(/a)') 'One-body parameters:'
738 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
739 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
743 C---------------------- GB or BP potential -----------------------------
744 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
745 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
747 C For the GB potential convert sigma'**2 into chi'
750 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
754 write (iout,'(/a/)') 'Parameters of the BP potential:'
755 write (iout,'(a/)') 'The epsilon array:'
756 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
757 write (iout,'(/a)') 'One-body parameters:'
758 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
760 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
761 & chip(i),alp(i),i=1,ntyp)
764 C--------------------- GBV potential -----------------------------------
765 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
766 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
767 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
769 write (iout,'(/a/)') 'Parameters of the GBV potential:'
770 write (iout,'(a/)') 'The epsilon array:'
771 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
772 write (iout,'(/a)') 'One-body parameters:'
773 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
774 & 's||/s_|_^2',' chip ',' alph '
775 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
776 & sigii(i),chip(i),alp(i),i=1,ntyp)
780 C-----------------------------------------------------------------------
781 C Calculate the "working" parameters of SC interactions.
789 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
790 sigma(j,i)=sigma(i,j)
791 rs0(i,j)=dwa16*sigma(i,j)
795 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
796 & 'Working parameters of the SC interactions:',
797 & ' a ',' b ',' augm ',' sigma ',' r0 ',
802 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
811 sigeps=dsign(1.0D0,epsij)
813 aa(i,j)=epsij*rrij*rrij
814 bb(i,j)=-sigeps*epsij*rrij
822 ratsig1=sigt2sq/sigt1sq
823 ratsig2=1.0D0/ratsig1
824 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
825 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
826 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
830 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
831 sigmaii(i,j)=rsum_max
832 sigmaii(j,i)=rsum_max
834 c sigmaii(i,j)=r0(i,j)
835 c sigmaii(j,i)=r0(i,j)
837 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
838 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
839 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
840 augm(i,j)=epsij*r_augm**(2*expon)
841 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
848 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
849 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
850 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
855 C Define the SC-p interaction constants
865 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
867 c aad(i,1)=0.3D0*4.0D0**12
868 C Following line for constants currently implemented
869 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
870 aad(i,1)=1.5D0*4.0D0**12
871 c aad(i,1)=0.17D0*5.6D0**12
873 C "Soft" SC-p repulsion
875 C Following line for constants currently implemented
876 c aad(i,1)=0.3D0*4.0D0**6
877 C "Hard" SC-p repulsion
878 bad(i,1)=3.0D0*4.0D0**6
879 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
888 C 8/9/01 Read the SC-p interaction constants from file
891 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
894 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
895 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
896 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
897 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
901 write (iout,*) "Parameters of SC-p interactions:"
903 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
904 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
909 C Define the constants of the disulfide bridge
913 c Old arbitrary potential - commented out.
918 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
919 c energy surface of diethyl disulfide.
920 c A. Liwo and U. Kozlowska, 11/24/03
930 write (iout,'(/a)') "Disulfide bridge parameters:"
931 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
932 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
933 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
934 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,