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 cc maxinter is maximum interaction sites
581 &nterm_sccor(i,j),nlor_sccor(i,j)
585 do k=1,nterm_sccor(i,j)
586 read (isccor,*) kk,v1sccor(k,l,i,j)
588 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
591 do k=1,nlor_sccor(i,j)
592 read (isccor,*) kk,vlor1sccor(k,i,j),
593 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
594 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
595 &(1+vlor3sccor(k,i,j)**2)
597 v0sccor(i,j)=v0ijsccor
604 write (iout,'(/a/)') 'Torsional constants:'
607 write (iout,*) 'ityp',i,' jtyp',j
608 write (iout,*) 'Fourier constants'
609 do k=1,nterm_sccor(i,j)
610 write (iout,'(2(1pe15.5))')
611 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
613 write (iout,*) 'Lorenz constants'
614 do k=1,nlor_sccor(i,j)
615 write (iout,'(3(1pe15.5))')
616 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
623 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
624 C interaction energy of the Gly, Ala, and Pro prototypes.
626 read (ifourier,*) nloctyp
629 read (ifourier,*) (b(ii,i),ii=1,13)
631 write (iout,*) 'Type',i
632 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
638 B1tilde(1,i) = b(3,i)
639 B1tilde(2,i) =-b(5,i)
640 B1tilde(1,-i) =-b(3,i)
641 B1tilde(2,-i) =b(5,i)
656 Ctilde(2,1,i)=-b(9,i)
658 Ctilde(1,1,-i)=b(7,i)
659 Ctilde(1,2,-i)=-b(9,i)
660 Ctilde(2,1,-i)=b(9,i)
661 Ctilde(2,2,-i)=b(7,i)
672 Dtilde(2,1,i)=-b(8,i)
674 Dtilde(1,1,-i)=b(6,i)
675 Dtilde(1,2,-i)=-b(8,i)
676 Dtilde(2,1,-i)=b(8,i)
677 Dtilde(2,2,-i)=b(6,i)
678 EE(1,1,i)= b(10,i)+b(11,i)
679 EE(2,2,i)=-b(10,i)+b(11,i)
680 EE(2,1,i)= b(12,i)-b(13,i)
681 EE(1,2,i)= b(12,i)+b(13,i)
682 EE(1,1,-i)= b(10,i)+b(11,i)
683 EE(2,2,-i)=-b(10,i)+b(11,i)
684 EE(2,1,-i)=-b(12,i)+b(13,i)
685 EE(1,2,-i)=-b(12,i)-b(13,i)
689 write (iout,*) 'Type',i
691 c write (iout,'(f10.5)') B1(:,i)
692 write(iout,*) B1(1,i),B1(2,i)
694 c write (iout,'(f10.5)') B2(:,i)
695 write(iout,*) B2(1,i),B2(2,i)
698 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
702 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
706 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
711 C Read electrostatic-interaction parameters
714 write (iout,'(/a)') 'Electrostatic interaction constants:'
715 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
716 & 'IT','JT','APP','BPP','AEL6','AEL3'
718 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
719 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
720 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
721 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
726 app (i,j)=epp(i,j)*rri*rri
727 bpp (i,j)=-2.0D0*epp(i,j)*rri
728 ael6(i,j)=elpp6(i,j)*4.2D0**6
729 ael3(i,j)=elpp3(i,j)*4.2D0**3
730 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
731 & ael6(i,j),ael3(i,j)
735 C Read side-chain interaction parameters.
737 read (isidep,*) ipot,expon
738 if (ipot.lt.1 .or. ipot.gt.5) then
739 write (iout,'(2a)') 'Error while reading SC interaction',
740 & 'potential file - unknown potential type.'
744 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
745 & ', exponents are ',expon,2*expon
746 goto (10,20,30,30,40) ipot
747 C----------------------- LJ potential ---------------------------------
748 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
750 write (iout,'(/a/)') 'Parameters of the LJ potential:'
751 write (iout,'(a/)') 'The epsilon array:'
752 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
753 write (iout,'(/a)') 'One-body parameters:'
754 write (iout,'(a,4x,a)') 'residue','sigma'
755 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
758 C----------------------- LJK potential --------------------------------
759 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
760 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
762 write (iout,'(/a/)') 'Parameters of the LJK potential:'
763 write (iout,'(a/)') 'The epsilon array:'
764 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
765 write (iout,'(/a)') 'One-body parameters:'
766 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
767 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
771 C---------------------- GB or BP potential -----------------------------
772 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
773 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
775 C For the GB potential convert sigma'**2 into chi'
778 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
782 write (iout,'(/a/)') 'Parameters of the BP potential:'
783 write (iout,'(a/)') 'The epsilon array:'
784 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
785 write (iout,'(/a)') 'One-body parameters:'
786 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
788 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
789 & chip(i),alp(i),i=1,ntyp)
792 C--------------------- GBV potential -----------------------------------
793 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
794 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
795 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
797 write (iout,'(/a/)') 'Parameters of the GBV 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,5a)') 'residue',' sigma ',' r0 ',
802 & 's||/s_|_^2',' chip ',' alph '
803 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
804 & sigii(i),chip(i),alp(i),i=1,ntyp)
808 C-----------------------------------------------------------------------
809 C Calculate the "working" parameters of SC interactions.
817 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
818 sigma(j,i)=sigma(i,j)
819 rs0(i,j)=dwa16*sigma(i,j)
823 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
824 & 'Working parameters of the SC interactions:',
825 & ' a ',' b ',' augm ',' sigma ',' r0 ',
830 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
839 sigeps=dsign(1.0D0,epsij)
841 aa(i,j)=epsij*rrij*rrij
842 bb(i,j)=-sigeps*epsij*rrij
850 ratsig1=sigt2sq/sigt1sq
851 ratsig2=1.0D0/ratsig1
852 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
853 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
854 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
858 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
859 sigmaii(i,j)=rsum_max
860 sigmaii(j,i)=rsum_max
862 c sigmaii(i,j)=r0(i,j)
863 c sigmaii(j,i)=r0(i,j)
865 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
866 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
867 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
868 augm(i,j)=epsij*r_augm**(2*expon)
869 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
876 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
877 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
878 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
883 C Define the SC-p interaction constants
893 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
895 c aad(i,1)=0.3D0*4.0D0**12
896 C Following line for constants currently implemented
897 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
898 aad(i,1)=1.5D0*4.0D0**12
899 c aad(i,1)=0.17D0*5.6D0**12
901 C "Soft" SC-p repulsion
903 C Following line for constants currently implemented
904 c aad(i,1)=0.3D0*4.0D0**6
905 C "Hard" SC-p repulsion
906 bad(i,1)=3.0D0*4.0D0**6
907 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
916 C 8/9/01 Read the SC-p interaction constants from file
919 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
922 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
923 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
924 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
925 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
929 write (iout,*) "Parameters of SC-p interactions:"
931 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
932 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
937 C Define the constants of the disulfide bridge
941 c Old arbitrary potential - commented out.
946 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
947 c energy surface of diethyl disulfide.
948 c A. Liwo and U. Kozlowska, 11/24/03
958 write (iout,'(/a)') "Disulfide bridge parameters:"
959 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
960 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
961 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
962 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,