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"/
22 character*1 toronelet(-2:2) /"p","a","G","A","P"/
24 dimension blower(3,3,maxlob)
25 double precision ip,mp
29 C Set LPRINT=.TRUE. for debugging
30 dwa16=2.0d0**(1.0d0/6.0d0)
33 C Assign virtual-bond length
38 c Read the virtual-bond parameters, masses, and moments of inertia
39 c and Stokes' radii of the peptide group and side chains
42 read (ibond,*) vbldp0,akp
45 read (ibond,*) vbldsc0(1,i),aksc(1,i)
50 dsc_inv(i)=1.0D0/dsc(i)
54 read (ibond,*) ijunk,vbldp0,akp,rjunk
56 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
62 dsc_inv(i)=1.0D0/dsc(i)
67 write(iout,'(/a/)')"Force constants virtual bonds:"
68 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
70 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
72 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
73 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
75 write (iout,'(13x,3f10.5)')
76 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
82 C Read the parameters of the probability distribution/energy expression
83 C of the virtual-bond valence angles theta
86 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
87 & (bthet(j,i,1,1),j=1,2)
88 read (ithep,*) (polthet(j,i),j=0,3)
89 read (ithep,*) (gthet(j,i),j=1,3)
90 read (ithep,*) theta0(i),sig0(i),sigc0(i)
94 athet(1,i,1,-1)=athet(1,i,1,1)
95 athet(2,i,1,-1)=athet(2,i,1,1)
96 bthet(1,i,1,-1)=-bthet(1,i,1,1)
97 bthet(2,i,1,-1)=-bthet(2,i,1,1)
98 athet(1,i,-1,1)=-athet(1,i,1,1)
99 athet(2,i,-1,1)=-athet(2,i,1,1)
100 bthet(1,i,-1,1)=bthet(1,i,1,1)
101 bthet(2,i,-1,1)=bthet(2,i,1,1)
105 athet(1,i,-1,-1)=athet(1,-i,1,1)
106 athet(2,i,-1,-1)=-athet(2,-i,1,1)
107 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
108 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
109 athet(1,i,-1,1)=athet(1,-i,1,1)
110 athet(2,i,-1,1)=-athet(2,-i,1,1)
111 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
112 bthet(2,i,-1,1)=bthet(2,-i,1,1)
113 athet(1,i,1,-1)=-athet(1,-i,1,1)
114 athet(2,i,1,-1)=athet(2,-i,1,1)
115 bthet(1,i,1,-1)=bthet(1,-i,1,1)
116 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
121 polthet(j,i)=polthet(j,-i)
124 gthet(j,i)=gthet(j,-i)
130 c & 'Parameters of the virtual-bond valence angles:'
131 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
132 c & ' ATHETA0 ',' A1 ',' A2 ',
135 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
136 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
138 c write (iout,'(/a/9x,5a/79(1h-))')
139 c & 'Parameters of the expression for sigma(theta_c):',
140 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
141 c & ' ALPH3 ',' SIGMA0C '
143 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
144 c & (polthet(j,i),j=0,3),sigc0(i)
146 c write (iout,'(/a/9x,5a/79(1h-))')
147 c & 'Parameters of the second gaussian:',
148 c & ' THETA0 ',' SIGMA0 ',' G1 ',
151 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
152 c & sig0(i),(gthet(j,i),j=1,3)
155 & 'Parameters of the virtual-bond valence angles:'
156 write (iout,'(/a/9x,5a/79(1h-))')
157 & 'Coefficients of expansion',
158 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
159 & ' b1*10^1 ',' b2*10^1 '
161 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
162 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
163 & (10*bthet(j,i,1,1),j=1,2)
165 write (iout,'(/a/9x,5a/79(1h-))')
166 & 'Parameters of the expression for sigma(theta_c):',
167 & ' alpha0 ',' alph1 ',' alph2 ',
168 & ' alhp3 ',' sigma0c '
170 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
171 & (polthet(j,i),j=0,3),sigc0(i)
173 write (iout,'(/a/9x,5a/79(1h-))')
174 & 'Parameters of the second gaussian:',
175 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
178 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
179 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
184 C Read the parameters of Utheta determined from ab initio surfaces
185 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
187 c write (iout,*) "tu dochodze"
188 read (ithep,*) nthetyp,ntheterm,ntheterm2,
189 & ntheterm3,nsingle,ndouble
190 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
191 read (ithep,*) (ithetyp(i),i=1,ntyp1)
193 ithetyp(i)=-ithetyp(-i)
195 c write (iout,*) "tu dochodze"
197 do i=-maxthetyp,maxthetyp
198 do j=-maxthetyp,maxthetyp
199 do k=-maxthetyp,maxthetyp
200 aa0thet(i,j,k,iblock)=0.0d0
202 aathet(l,i,j,k,iblock)=0.0d0
206 bbthet(m,l,i,j,k,iblock)=0.0d0
207 ccthet(m,l,i,j,k,iblock)=0.0d0
208 ddthet(m,l,i,j,k,iblock)=0.0d0
209 eethet(m,l,i,j,k,iblock)=0.0d0
215 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
216 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
226 do j=-nthetyp,nthetyp
227 do k=-nthetyp,nthetyp
228 read (ithep,'(6a)') res1
229 read (ithep,*) aa0thet(i,j,k,iblock)
230 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
232 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
233 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
234 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
235 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
238 & (((ffthet(llll,lll,ll,i,j,k,iblock),
239 & ffthet(lll,llll,ll,i,j,k,iblock),
240 & ggthet(llll,lll,ll,i,j,k,iblock)
241 & ,ggthet(lll,llll,ll,i,j,k,iblock),
242 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
247 C For dummy ends assign glycine-type coefficients of theta-only terms; the
248 C coefficients of theta-and-gamma-dependent terms are zero.
253 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
254 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
256 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
257 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
260 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
262 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
265 C Substitution for D aminoacids from symmetry.
268 do j=-nthetyp,nthetyp
269 do k=-nthetyp,nthetyp
270 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
272 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
276 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
277 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
278 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
279 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
285 ffthet(llll,lll,ll,i,j,k,iblock)=
286 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
287 ffthet(lll,llll,ll,i,j,k,iblock)=
288 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
289 ggthet(llll,lll,ll,i,j,k,iblock)=
290 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
291 ggthet(lll,llll,ll,i,j,k,iblock)=
292 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
302 C Control printout of the coefficients of virtual-bond-angle potentials
305 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
310 write (iout,'(//4a)')
311 & 'Type ',onelett(i),onelett(j),onelett(k)
312 write (iout,'(//a,10x,a)') " l","a[l]"
313 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
314 write (iout,'(i2,1pe15.5)')
315 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
317 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
318 & "b",l,"c",l,"d",l,"e",l
320 write (iout,'(i2,4(1pe15.5))') m,
321 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
322 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
326 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
327 & "f+",l,"f-",l,"g+",l,"g-",l
330 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
331 & ffthet(n,m,l,i,j,k,iblock),
332 & ffthet(m,n,l,i,j,k,iblock),
333 & ggthet(n,m,l,i,j,k,iblock),
334 & ggthet(m,n,l,i,j,k,iblock)
348 C Read the parameters of the probability distribution/energy expression
349 C of the side chains.
352 cc write (iout,*) "tu dochodze",i
353 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
357 dsc_inv(i)=1.0D0/dsc(i)
368 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
369 censc(1,1,-i)=censc(1,1,i)
370 censc(2,1,-i)=censc(2,1,i)
371 censc(3,1,-i)=-censc(3,1,i)
373 read (irotam,*) bsc(j,i)
374 read (irotam,*) (censc(k,j,i),k=1,3),
375 & ((blower(k,l,j),l=1,k),k=1,3)
376 censc(1,j,-i)=censc(1,j,i)
377 censc(2,j,-i)=censc(2,j,i)
378 censc(3,j,-i)=-censc(3,j,i)
379 C BSC is amplitude of Gaussian
386 akl=akl+blower(k,m,j)*blower(l,m,j)
390 if (((k.eq.3).and.(l.ne.3))
391 & .or.((l.eq.3).and.(k.ne.3))) then
392 gaussc(k,l,j,-i)=-akl
393 gaussc(l,k,j,-i)=-akl
405 write (iout,'(/a)') 'Parameters of side-chain local geometry'
409 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
410 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
411 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
412 c write (iout,'(a,f10.4,4(16x,f10.4))')
413 c & 'Center ',(bsc(j,i),j=1,nlobi)
414 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
415 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
416 & 'log h',(bsc(j,i),j=1,nlobi)
417 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
418 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
425 c blower(k,l,j)=gaussc(ind,j,i)
430 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
431 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
438 C Read scrot parameters for potentials determined from all-atom AM1 calculations
439 C added by Urszula Kozlowska 07/11/2007
447 read(irotam,*) sc_parmin(j,i)
455 C Read torsional parameters in old format
457 read (itorp,*) ntortyp,nterm_old
458 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
459 read (itorp,*) (itortyp(i),i=1,ntyp)
464 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
470 write (iout,'(/a/)') 'Torsional constants:'
473 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
474 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
482 C Read torsional parameters
484 read (itorp,*) ntortyp
485 read (itorp,*) (itortyp(i),i=1,ntyp)
486 write (iout,*) 'ntortyp',ntortyp
489 itortyp(i)=-itortyp(-i)
491 c write (iout,*) 'ntortyp',ntortyp
493 do j=-ntortyp+1,ntortyp-1
494 read (itorp,*) nterm(i,j,iblock),
496 nterm(-i,-j,iblock)=nterm(i,j,iblock)
497 nlor(-i,-j,iblock)=nlor(i,j,iblock)
500 do k=1,nterm(i,j,iblock)
501 read (itorp,*) kk,v1(k,i,j,iblock),
503 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
504 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
505 v0ij=v0ij+si*v1(k,i,j,iblock)
508 do k=1,nlor(i,j,iblock)
509 read (itorp,*) kk,vlor1(k,i,j),
510 & vlor2(k,i,j),vlor3(k,i,j)
511 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
514 v0(-i,-j,iblock)=v0ij
520 write (iout,'(/a/)') 'Torsional constants:'
523 write (iout,*) 'ityp',i,' jtyp',j
524 write (iout,*) 'Fourier constants'
526 do k=1,nterm(i,j,iblock)
527 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
530 write (iout,*) 'Lorenz constants'
531 do k=1,nlor(i,j,iblock)
532 write (iout,'(3(1pe15.5))')
533 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
540 C 6/23/01 Read parameters for double torsionals
544 do j=-ntortyp+1,ntortyp-1
545 do k=-ntortyp+1,ntortyp-1
546 read (itordp,'(3a1)') t1,t2,t3
547 c write (iout,*) "OK onelett",
550 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
551 & .or. t3.ne.toronelet(k)) then
552 write (iout,*) "Error in double torsional parameter file",
555 call MPI_Finalize(Ierror)
557 stop "Error in double torsional parameter file"
559 read (itordp,*) ntermd_1(i,j,k,iblock),
560 & ntermd_2(i,j,k,iblock)
561 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
562 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
563 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
564 & ntermd_1(i,j,k,iblock))
565 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
566 & ntermd_1(i,j,k,iblock))
567 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
568 & ntermd_1(i,j,k,iblock))
569 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
570 & ntermd_1(i,j,k,iblock))
571 C Martix of D parameters for one dimesional foureir series
572 do l=1,ntermd_1(i,j,k,iblock)
573 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
574 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
575 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
576 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
577 c write(iout,*) "whcodze" ,
578 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
580 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
581 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
582 & v2s(m,l,i,j,k,iblock),
583 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
584 C Martix of D parameters for two dimesional fourier series
585 do l=1,ntermd_2(i,j,k,iblock)
587 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
588 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
589 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
590 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
599 write (iout,*) 'Constants for double torsionals'
602 do j=-ntortyp+1,ntortyp-1
603 do k=-ntortyp+1,ntortyp-1
604 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
605 & ' nsingle',ntermd_1(i,j,k,iblock),
606 & ' ndouble',ntermd_2(i,j,k,iblock)
608 write (iout,*) 'Single angles:'
609 do l=1,ntermd_1(i,j,k,iblock)
610 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
611 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
612 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
613 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
616 write (iout,*) 'Pairs of angles:'
617 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
618 do l=1,ntermd_2(i,j,k,iblock)
619 write (iout,'(i5,20f10.5)')
620 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
623 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
624 do l=1,ntermd_2(i,j,k,iblock)
625 write (iout,'(i5,20f10.5)')
626 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
627 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
636 C Read of Side-chain backbone correlation parameters
637 C Modified 11 May 2012 by Adasko
640 read (isccor,*) nsccortyp
641 read (isccor,*) (isccortyp(i),i=1,ntyp)
643 isccortyp(i)=-isccortyp(-i)
645 iscprol=isccortyp(20)
646 c write (iout,*) 'ntortyp',ntortyp
648 cc maxinter is maximum interaction sites
653 &nterm_sccor(i,j),nlor_sccor(i,j)
654 write (iout,*) nterm_sccor(i,j)
660 nterm_sccor(-i,j)=nterm_sccor(i,j)
661 nterm_sccor(-i,-j)=nterm_sccor(i,j)
662 nterm_sccor(i,-j)=nterm_sccor(i,j)
663 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
664 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
665 do k=1,nterm_sccor(i,j)
666 read (isccor,*) kk,v1sccor(k,l,i,j)
668 if (j.eq.iscprol) then
669 if (i.eq.isccortyp(10)) then
670 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
671 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
673 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
674 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
675 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
676 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
677 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
678 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
679 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
680 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
683 if (i.eq.isccortyp(10)) then
684 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
685 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
687 if (j.eq.isccortyp(10)) then
688 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
689 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
691 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
692 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
693 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
694 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
695 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
696 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
700 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
701 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
702 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
703 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
706 do k=1,nlor_sccor(i,j)
707 read (isccor,*) kk,vlor1sccor(k,i,j),
708 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
709 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
710 &(1+vlor3sccor(k,i,j)**2)
712 v0sccor(l,i,j)=v0ijsccor
713 v0sccor(l,-i,j)=v0ijsccor1
714 v0sccor(l,i,-j)=v0ijsccor2
715 v0sccor(l,-i,-j)=v0ijsccor3
721 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
724 write (iout,*) 'ityp',i,' jtyp',j
725 write (iout,*) 'Fourier constants'
727 do k=1,nterm_sccor(i,j)
728 write (iout,'(2(1pe15.5))')
729 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
732 write (iout,*) 'Lorenz constants'
733 do k=1,nlor_sccor(i,j)
734 write (iout,'(3(1pe15.5))')
735 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
741 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
742 C interaction energy of the Gly, Ala, and Pro prototypes.
744 write (iout,*) "Reading Fourier"
746 read (ifourier,*) nloctyp
749 read (ifourier,*) (b(ii,i),ii=1,13)
751 write (iout,*) 'Type',i
752 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
754 write (iout,*) "Reading Fourier 1"
757 read (ifourier,*) (bnew1(ii,1,i),ii=1,3)
758 read (ifourier,*) (bnew2(ii,1,i),ii=1,3)
759 read (ifourier,*) (bnew1(ii,2,i),ii=1,1)
760 read (ifourier,*) (bnew2(ii,2,i),ii=1,1)
761 read (ifourier,*) (eenew(ii,i),ii=1,1)
772 Ctilde(2,1,i)=-b(9,i)
774 Ctilde(1,1,-i)=b(7,i)
775 Ctilde(1,2,-i)=-b(9,i)
776 Ctilde(2,1,-i)=b(9,i)
777 Ctilde(2,2,-i)=b(7,i)
789 Dtilde(2,1,i)=-b(8,i)
791 Dtilde(1,1,-i)=b(6,i)
792 Dtilde(1,2,-i)=-b(8,i)
793 Dtilde(2,1,-i)=b(8,i)
794 Dtilde(2,2,-i)=b(6,i)
796 EEold(1,1,i)= b(10,i)+b(11,i)
797 eeold(1,1,i)= eenew(1,i)
798 EEold(2,2,i)=-b(10,i)+b(11,i)
799 EEold(2,1,i)= b(12,i)-b(13,i)
800 EEold(1,2,i)= b(12,i)+b(13,i)
801 EEold(1,1,-i)= b(10,i)+b(11,i)
802 EEold(2,2,-i)=-b(10,i)+b(11,i)
803 EEold(2,1,-i)=-b(12,i)+b(13,i)
804 EEold(1,2,-i)=-b(12,i)-b(13,i)
805 c write(iout,*) "TU DOCHODZE"
810 write (iout,*) 'Type',i
811 write (iout,*) 'BNEW1(1)'
812 write (iout,*) (BNEW1(ii,1,i),ii=1,3)
813 write (iout,*) 'BNEW1(2)'
814 write (iout,*) BNEW1(1,2,i)
815 write (iout,*) 'BNEW2(1)'
816 write (iout,*) (BNEW2(ii,1,i),ii=1,3)
817 write (iout,*) 'BNEW2(2)'
818 write (iout,*) BNEW2(1,2,i)
821 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
825 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
829 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
840 B1tilde(1,i) = b(3,i)
841 B1tilde(2,i) =-b(5,i)
842 B1tilde(1,-i) =-b(3,i)
843 B1tilde(2,-i) =b(5,i)
867 Ctilde(2,1,i)=-b(9,i)
869 Ctilde(1,1,-i)=b(7,i)
870 Ctilde(1,2,-i)=-b(9,i)
871 Ctilde(2,1,-i)=b(9,i)
872 Ctilde(2,2,-i)=b(7,i)
874 c Ctilde(1,1,i)=0.0d0
875 c Ctilde(1,2,i)=0.0d0
876 c Ctilde(2,1,i)=0.0d0
877 c Ctilde(2,2,i)=0.0d0
892 Dtilde(2,1,i)=-b(8,i)
894 Dtilde(1,1,-i)=b(6,i)
895 Dtilde(1,2,-i)=-b(8,i)
896 Dtilde(2,1,-i)=b(8,i)
897 Dtilde(2,2,-i)=b(6,i)
899 c Dtilde(1,1,i)=0.0d0
900 c Dtilde(1,2,i)=0.0d0
901 c Dtilde(2,1,i)=0.0d0
902 c Dtilde(2,2,i)=0.0d0
903 EE(1,1,i)= b(10,i)+b(11,i)
904 EE(2,2,i)=-b(10,i)+b(11,i)
905 EE(2,1,i)= b(12,i)-b(13,i)
906 EE(1,2,i)= b(12,i)+b(13,i)
907 EE(1,1,-i)= b(10,i)+b(11,i)
908 EE(2,2,-i)=-b(10,i)+b(11,i)
909 EE(2,1,-i)=-b(12,i)+b(13,i)
910 EE(1,2,-i)=-b(12,i)-b(13,i)
916 c ee(2,1,i)=ee(1,2,i)
921 write (iout,*) 'Type',i
923 c write (iout,'(f10.5)') B1(:,i)
924 write(iout,*) B1(1,i),B1(2,i)
926 c write (iout,'(f10.5)') B2(:,i)
927 write(iout,*) B2(1,i),B2(2,i)
930 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
934 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
938 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
944 C Read electrostatic-interaction parameters
947 write (iout,'(/a)') 'Electrostatic interaction constants:'
948 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
949 & 'IT','JT','APP','BPP','AEL6','AEL3'
951 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
952 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
953 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
954 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
959 app (i,j)=epp(i,j)*rri*rri
960 bpp (i,j)=-2.0D0*epp(i,j)*rri
961 ael6(i,j)=elpp6(i,j)*4.2D0**6
962 ael3(i,j)=elpp3(i,j)*4.2D0**3
964 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
965 & ael6(i,j),ael3(i,j)
970 C Read side-chain interaction parameters.
972 read (isidep,*) ipot,expon
973 if (ipot.lt.1 .or. ipot.gt.5) then
974 write (iout,'(2a)') 'Error while reading SC interaction',
975 & 'potential file - unknown potential type.'
979 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
980 & ', exponents are ',expon,2*expon
981 goto (10,20,30,30,40) ipot
982 C----------------------- LJ potential ---------------------------------
983 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
985 write (iout,'(/a/)') 'Parameters of the LJ potential:'
986 write (iout,'(a/)') 'The epsilon array:'
987 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
988 write (iout,'(/a)') 'One-body parameters:'
989 write (iout,'(a,4x,a)') 'residue','sigma'
990 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
993 C----------------------- LJK potential --------------------------------
994 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
995 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
997 write (iout,'(/a/)') 'Parameters of the LJK potential:'
998 write (iout,'(a/)') 'The epsilon array:'
999 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1000 write (iout,'(/a)') 'One-body parameters:'
1001 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1002 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1006 C---------------------- GB or BP potential -----------------------------
1008 read (isidep,*)(eps(i,j),j=i,ntyp)
1010 read (isidep,*)(sigma0(i),i=1,ntyp)
1011 read (isidep,*)(sigii(i),i=1,ntyp)
1012 read (isidep,*)(chip(i),i=1,ntyp)
1013 read (isidep,*)(alp(i),i=1,ntyp)
1014 C For the GB potential convert sigma'**2 into chi'
1017 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1021 write (iout,'(/a/)') 'Parameters of the BP potential:'
1022 write (iout,'(a/)') 'The epsilon array:'
1023 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1024 write (iout,'(/a)') 'One-body parameters:'
1025 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1027 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1028 & chip(i),alp(i),i=1,ntyp)
1031 C--------------------- GBV potential -----------------------------------
1032 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1033 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1034 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1036 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1037 write (iout,'(a/)') 'The epsilon array:'
1038 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1039 write (iout,'(/a)') 'One-body parameters:'
1040 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1041 & 's||/s_|_^2',' chip ',' alph '
1042 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1043 & sigii(i),chip(i),alp(i),i=1,ntyp)
1047 C-----------------------------------------------------------------------
1048 C Calculate the "working" parameters of SC interactions.
1056 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1057 sigma(j,i)=sigma(i,j)
1058 rs0(i,j)=dwa16*sigma(i,j)
1062 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1063 & 'Working parameters of the SC interactions:',
1064 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1069 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1078 sigeps=dsign(1.0D0,epsij)
1080 aa(i,j)=epsij*rrij*rrij
1081 bb(i,j)=-sigeps*epsij*rrij
1085 sigt1sq=sigma0(i)**2
1086 sigt2sq=sigma0(j)**2
1089 ratsig1=sigt2sq/sigt1sq
1090 ratsig2=1.0D0/ratsig1
1091 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1092 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1093 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1097 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1098 sigmaii(i,j)=rsum_max
1099 sigmaii(j,i)=rsum_max
1101 c sigmaii(i,j)=r0(i,j)
1102 c sigmaii(j,i)=r0(i,j)
1104 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1105 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1106 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1107 augm(i,j)=epsij*r_augm**(2*expon)
1108 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1115 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1116 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1117 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1122 C Define the SC-p interaction constants
1126 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1128 c aad(i,1)=0.3D0*4.0D0**12
1129 C Following line for constants currently implemented
1130 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1131 aad(i,1)=1.5D0*4.0D0**12
1132 c aad(i,1)=0.17D0*5.6D0**12
1134 C "Soft" SC-p repulsion
1136 C Following line for constants currently implemented
1137 c aad(i,1)=0.3D0*4.0D0**6
1138 C "Hard" SC-p repulsion
1139 bad(i,1)=3.0D0*4.0D0**6
1140 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1149 C 8/9/01 Read the SC-p interaction constants from file
1152 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1155 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1156 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1157 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1158 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1162 write (iout,*) "Parameters of SC-p interactions:"
1164 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1165 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1170 C Define the constants of the disulfide bridge
1174 c Old arbitrary potential - commented out.
1179 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1180 c energy surface of diethyl disulfide.
1181 c A. Liwo and U. Kozlowska, 11/24/03
1192 write (iout,'(/a)') "Disulfide bridge parameters:"
1193 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1194 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1195 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1196 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,