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'
20 include 'COMMON.SHIELD'
21 include 'COMMON.CONTROL'
24 character*1 onelett(4) /"G","A","P","D"/
25 character*1 toronelet(-2:2) /"p","a","G","A","P"/
27 dimension blower(3,3,maxlob)
28 double precision ip,mp
29 character*3 lancuch,ucase
33 write (iout,*) "PARMREAD tor_mode",tor_mode
34 call getenv("PRINT_PARM",lancuch)
35 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
36 write (iout,*) "lprint ",lprint
37 C Set LPRINT=.TRUE. for debugging
38 dwa16=2.0d0**(1.0d0/6.0d0)
40 C Assign virtual-bond length
45 read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp
48 read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i)
53 dsc_inv(i)=1.0D0/dsc(i)
57 read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk
59 read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
60 & aksc(j,i),abond0(j,i),j=1,nbondterm(i))
65 dsc_inv(i)=1.0D0/dsc(i)
70 write(iout,'(/a/)')"Force constants virtual bonds:"
71 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
73 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
75 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
76 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
78 write (iout,'(13x,3f10.5)')
79 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
83 read(iliptranpar,*,end=1161,err=1161) pepliptran
85 read(iliptranpar,*,end=1161,err=1161) liptranene(i)
90 C Read the parameters of the probability distribution/energy expression
91 C of the virtual-bond valence angles theta
94 read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
95 & (bthet(j,i,1,1),j=1,2)
96 read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
97 read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
98 read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
102 athet(1,i,1,-1)=athet(1,i,1,1)
103 athet(2,i,1,-1)=athet(2,i,1,1)
104 bthet(1,i,1,-1)=-bthet(1,i,1,1)
105 bthet(2,i,1,-1)=-bthet(2,i,1,1)
106 athet(1,i,-1,1)=-athet(1,i,1,1)
107 athet(2,i,-1,1)=-athet(2,i,1,1)
108 bthet(1,i,-1,1)=bthet(1,i,1,1)
109 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)
117 athet(1,i,-1,1)=athet(1,-i,1,1)
118 athet(2,i,-1,1)=-athet(2,-i,1,1)
119 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
120 bthet(2,i,-1,1)=bthet(2,-i,1,1)
121 athet(1,i,1,-1)=-athet(1,-i,1,1)
122 athet(2,i,1,-1)=athet(2,-i,1,1)
123 bthet(1,i,1,-1)=bthet(1,-i,1,1)
124 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
129 polthet(j,i)=polthet(j,-i)
132 gthet(j,i)=gthet(j,-i)
138 c & 'Parameters of the virtual-bond valence angles:'
139 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
140 c & ' ATHETA0 ',' A1 ',' A2 ',
143 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
144 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
146 c write (iout,'(/a/9x,5a/79(1h-))')
147 c & 'Parameters of the expression for sigma(theta_c):',
148 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
149 c & ' ALPH3 ',' SIGMA0C '
151 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
152 c & (polthet(j,i),j=0,3),sigc0(i)
154 c write (iout,'(/a/9x,5a/79(1h-))')
155 c & 'Parameters of the second gaussian:',
156 c & ' THETA0 ',' SIGMA0 ',' G1 ',
159 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
160 c & sig0(i),(gthet(j,i),j=1,3)
163 & 'Parameters of the virtual-bond valence angles:'
164 write (iout,'(/a/9x,5a/79(1h-))')
165 & 'Coefficients of expansion',
166 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
167 & ' b1*10^1 ',' b2*10^1 '
169 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
170 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
171 & (10*bthet(j,i,1,1),j=1,2)
173 write (iout,'(/a/9x,5a/79(1h-))')
174 & 'Parameters of the expression for sigma(theta_c):',
175 & ' alpha0 ',' alph1 ',' alph2 ',
176 & ' alhp3 ',' sigma0c '
178 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
179 & (polthet(j,i),j=0,3),sigc0(i)
181 write (iout,'(/a/9x,5a/79(1h-))')
182 & 'Parameters of the second gaussian:',
183 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
186 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
187 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
191 IF (tor_mode.eq.0) THEN
193 C Read the parameters of Utheta determined from ab initio surfaces
194 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
196 read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
197 & ntheterm3,nsingle,ndouble
198 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
199 read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
201 ithetyp(i)=-ithetyp(-i)
204 do i=-maxthetyp,maxthetyp
205 do j=-maxthetyp,maxthetyp
206 do k=-maxthetyp,maxthetyp
207 aa0thet(i,j,k,iblock)=0.0d0
209 aathet(l,i,j,k,iblock)=0.0d0
213 bbthet(m,l,i,j,k,iblock)=0.0d0
214 ccthet(m,l,i,j,k,iblock)=0.0d0
215 ddthet(m,l,i,j,k,iblock)=0.0d0
216 eethet(m,l,i,j,k,iblock)=0.0d0
222 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
223 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
231 C write (iout,*) "KURWA1"
234 do j=-nthetyp,nthetyp
235 do k=-nthetyp,nthetyp
236 read (ithep,'(6a)',end=111,err=111) res1
237 write(iout,*) res1,i,j,k
238 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
239 read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
241 read (ithep,*,end=111,err=111)
242 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
243 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
244 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
245 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
247 read (ithep,*,end=111,err=111)
248 & (((ffthet(llll,lll,ll,i,j,k,iblock),
249 & ffthet(lll,llll,ll,i,j,k,iblock),
250 & ggthet(llll,lll,ll,i,j,k,iblock)
251 & ,ggthet(lll,llll,ll,i,j,k,iblock),
252 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
256 C write(iout,*) "KURWA1.1"
258 C For dummy ends assign glycine-type coefficients of theta-only terms; the
259 C coefficients of theta-and-gamma-dependent terms are zero.
264 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
265 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
267 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
268 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
271 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
273 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
276 C write(iout,*) "KURWA1.5"
277 C Substitution for D aminoacids from symmetry.
280 do j=-nthetyp,nthetyp
281 do k=-nthetyp,nthetyp
282 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
284 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
288 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
289 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
290 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
291 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
297 ffthet(llll,lll,ll,i,j,k,iblock)=
298 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
299 ffthet(lll,llll,ll,i,j,k,iblock)=
300 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
301 ggthet(llll,lll,ll,i,j,k,iblock)=
302 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
303 ggthet(lll,llll,ll,i,j,k,iblock)=
304 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
314 C Control printout of the coefficients of virtual-bond-angle potentials
317 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
322 write (iout,'(//4a)')
323 & 'Type ',onelett(i),onelett(j),onelett(k)
324 write (iout,'(//a,10x,a)') " l","a[l]"
325 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
326 write (iout,'(i2,1pe15.5)')
327 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
329 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
330 & "b",l,"c",l,"d",l,"e",l
332 write (iout,'(i2,4(1pe15.5))') m,
333 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
334 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
338 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
339 & "f+",l,"f-",l,"g+",l,"g-",l
342 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
343 & ffthet(n,m,l,i,j,k,iblock),
344 & ffthet(m,n,l,i,j,k,iblock),
345 & ggthet(n,m,l,i,j,k,iblock),
346 & ggthet(m,n,l,i,j,k,iblock)
359 C here will be the apropriate recalibrating for D-aminoacid
360 read (ithep,*,end=111,err=111) nthetyp
361 do i=-nthetyp+1,nthetyp-1
362 read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
363 do j=0,nbend_kcc_Tb(i)
364 read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
369 & "Parameters of the valence-only potentials"
370 do i=-nthetyp+1,nthetyp-1
371 write (iout,'(2a)') "Type ",toronelet(i)
372 do k=0,nbend_kcc_Tb(i)
373 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
382 C write(iout,*) 'KURWA2'
385 C Read the parameters of the probability distribution/energy expression
386 C of the side chains.
389 cc write (iout,*) "tu dochodze",i
390 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
394 dsc_inv(i)=1.0D0/dsc(i)
405 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
406 & ((blower(k,l,1),l=1,k),k=1,3)
407 censc(1,1,-i)=censc(1,1,i)
408 censc(2,1,-i)=censc(2,1,i)
409 censc(3,1,-i)=-censc(3,1,i)
411 read (irotam,*,end=112,err=112) bsc(j,i)
412 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
413 & ((blower(k,l,j),l=1,k),k=1,3)
414 censc(1,j,-i)=censc(1,j,i)
415 censc(2,j,-i)=censc(2,j,i)
416 censc(3,j,-i)=-censc(3,j,i)
417 C BSC is amplitude of Gaussian
424 akl=akl+blower(k,m,j)*blower(l,m,j)
428 if (((k.eq.3).and.(l.ne.3))
429 & .or.((l.eq.3).and.(k.ne.3))) then
430 gaussc(k,l,j,-i)=-akl
431 gaussc(l,k,j,-i)=-akl
443 write (iout,'(/a)') 'Parameters of side-chain local geometry'
447 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
448 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
449 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
450 c write (iout,'(a,f10.4,4(16x,f10.4))')
451 c & 'Center ',(bsc(j,i),j=1,nlobi)
452 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
453 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
454 & 'log h',(bsc(j,i),j=1,nlobi)
455 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
456 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
463 c blower(k,l,j)=gaussc(ind,j,i)
468 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
469 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
476 C Read scrot parameters for potentials determined from all-atom AM1 calculations
477 C added by Urszula Kozlowska 07/11/2007
480 read (irotam,*,end=112,err=112)
482 read (irotam,*,end=112,err=112)
485 read(irotam,*,end=112,err=112) sc_parmin(j,i)
492 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
493 C interaction energy of the Gly, Ala, and Pro prototypes.
495 read (ifourier,*,end=115,err=115) nloctyp
496 SPLIT_FOURIERTOR = nloctyp.lt.0
497 nloctyp = iabs(nloctyp)
499 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
500 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
501 itype2loc(ntyp1)=nloctyp
502 iloctyp(nloctyp)=ntyp1
504 itype2loc(-i)=-itype2loc(i)
513 iloctyp(-i)=-iloctyp(i)
515 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
516 c write (iout,*) "nloctyp",nloctyp,
517 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
520 c write (iout,*) "NEWCORR",i
521 read (ifourier,*,end=115,err=115)
524 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
527 c write (iout,*) "NEWCORR BNEW1"
528 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
531 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
534 c write (iout,*) "NEWCORR BNEW2"
535 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
537 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
538 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
540 c write (iout,*) "NEWCORR CCNEW"
541 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
543 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
544 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
546 c write (iout,*) "NEWCORR DDNEW"
547 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
551 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
555 c write (iout,*) "NEWCORR EENEW1"
556 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
558 read (ifourier,*,end=115,err=115) e0new(ii,i)
560 c write (iout,*) (e0new(ii,i),ii=1,3)
562 c write (iout,*) "NEWCORR EENEW"
565 ccnew(ii,1,i)=ccnew(ii,1,i)/2
566 ccnew(ii,2,i)=ccnew(ii,2,i)/2
567 ddnew(ii,1,i)=ddnew(ii,1,i)/2
568 ddnew(ii,2,i)=ddnew(ii,2,i)/2
573 bnew1(ii,1,-i)= bnew1(ii,1,i)
574 bnew1(ii,2,-i)=-bnew1(ii,2,i)
575 bnew2(ii,1,-i)= bnew2(ii,1,i)
576 bnew2(ii,2,-i)=-bnew2(ii,2,i)
579 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
580 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
581 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
582 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
583 ccnew(ii,1,-i)=ccnew(ii,1,i)
584 ccnew(ii,2,-i)=-ccnew(ii,2,i)
585 ddnew(ii,1,-i)=ddnew(ii,1,i)
586 ddnew(ii,2,-i)=-ddnew(ii,2,i)
588 e0new(1,-i)= e0new(1,i)
589 e0new(2,-i)=-e0new(2,i)
590 e0new(3,-i)=-e0new(3,i)
592 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
593 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
594 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
595 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
599 write (iout,'(a)') "Coefficients of the multibody terms"
600 do i=-nloctyp+1,nloctyp-1
601 write (iout,*) "Type: ",onelet(iloctyp(i))
602 write (iout,*) "Coefficients of the expansion of B1"
604 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
606 write (iout,*) "Coefficients of the expansion of B2"
608 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
610 write (iout,*) "Coefficients of the expansion of C"
611 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
612 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
613 write (iout,*) "Coefficients of the expansion of D"
614 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
615 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
616 write (iout,*) "Coefficients of the expansion of E"
617 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
620 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
625 IF (SPLIT_FOURIERTOR) THEN
627 c write (iout,*) "NEWCORR TOR",i
628 read (ifourier,*,end=115,err=115)
631 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
634 c write (iout,*) "NEWCORR BNEW1 TOR"
635 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
638 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
641 c write (iout,*) "NEWCORR BNEW2 TOR"
642 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
644 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
645 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
647 c write (iout,*) "NEWCORR CCNEW TOR"
648 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
650 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
651 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
653 c write (iout,*) "NEWCORR DDNEW TOR"
654 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
658 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
662 c write (iout,*) "NEWCORR EENEW1 TOR"
663 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
665 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
667 c write (iout,*) (e0newtor(ii,i),ii=1,3)
669 c write (iout,*) "NEWCORR EENEW TOR"
672 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
673 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
674 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
675 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
680 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
681 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
682 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
683 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
686 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
687 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
688 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
689 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
691 e0newtor(1,-i)= e0newtor(1,i)
692 e0newtor(2,-i)=-e0newtor(2,i)
693 e0newtor(3,-i)=-e0newtor(3,i)
695 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
696 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
697 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
698 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
703 & "Single-body coefficients of the torsional potentials"
704 do i=-nloctyp+1,nloctyp-1
705 write (iout,*) "Type: ",onelet(iloctyp(i))
706 write (iout,*) "Coefficients of the expansion of B1tor"
708 write (iout,'(3hB1(,i1,1h),3f10.5)')
709 & j,(bnew1tor(k,j,i),k=1,3)
711 write (iout,*) "Coefficients of the expansion of B2tor"
713 write (iout,'(3hB2(,i1,1h),3f10.5)')
714 & j,(bnew2tor(k,j,i),k=1,3)
716 write (iout,*) "Coefficients of the expansion of Ctor"
717 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
718 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
719 write (iout,*) "Coefficients of the expansion of Dtor"
720 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
721 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
722 write (iout,*) "Coefficients of the expansion of Etor"
723 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
726 write (iout,'(1hE,2i1,2f10.5)')
727 & j,k,(eenewtor(l,j,k,i),l=1,2)
733 do i=-nloctyp+1,nloctyp-1
736 bnew1tor(ii,j,i)=bnew1(ii,j,i)
741 bnew2tor(ii,j,i)=bnew2(ii,j,i)
745 ccnewtor(ii,1,i)=ccnew(ii,1,i)
746 ccnewtor(ii,2,i)=ccnew(ii,2,i)
747 ddnewtor(ii,1,i)=ddnew(ii,1,i)
748 ddnewtor(ii,2,i)=ddnew(ii,2,i)
751 ENDIF !SPLIT_FOURIER_TOR
754 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
756 read (ifourier,*,end=115,err=115)
757 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
759 write (iout,*) 'Type ',onelet(iloctyp(i))
760 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
774 c B1tilde(1,i) = b(3)
775 c B1tilde(2,i) =-b(5)
776 c B1tilde(1,-i) =-b(3)
777 c B1tilde(2,-i) =b(5)
784 cc B1tilde(1,i) = b(3,i)
785 cc B1tilde(2,i) =-b(5,i)
786 C B1tilde(1,-i) =-b(3,i)
787 C B1tilde(2,-i) =b(5,i)
788 cc b1tilde(1,i)=0.0d0
789 cc b1tilde(2,i)=0.0d0
801 CCold(1,1,-i)= b(7,i)
802 CCold(2,2,-i)=-b(7,i)
803 CCold(2,1,-i)=-b(9,i)
804 CCold(1,2,-i)=-b(9,i)
809 c Ctilde(1,1,i)= CCold(1,1,i)
810 c Ctilde(1,2,i)= CCold(1,2,i)
811 c Ctilde(2,1,i)=-CCold(2,1,i)
812 c Ctilde(2,2,i)=-CCold(2,2,i)
814 c Ctilde(1,1,i)=0.0d0
815 c Ctilde(1,2,i)=0.0d0
816 c Ctilde(2,1,i)=0.0d0
817 c Ctilde(2,2,i)=0.0d0
822 DDold(1,1,-i)= b(6,i)
823 DDold(2,2,-i)=-b(6,i)
824 DDold(2,1,-i)=-b(8,i)
825 DDold(1,2,-i)=-b(8,i)
830 c Dtilde(1,1,i)= DD(1,1,i)
831 c Dtilde(1,2,i)= DD(1,2,i)
832 c Dtilde(2,1,i)=-DD(2,1,i)
833 c Dtilde(2,2,i)=-DD(2,2,i)
835 c Dtilde(1,1,i)=0.0d0
836 c Dtilde(1,2,i)=0.0d0
837 c Dtilde(2,1,i)=0.0d0
838 c Dtilde(2,2,i)=0.0d0
839 EEold(1,1,i)= b(10,i)+b(11,i)
840 EEold(2,2,i)=-b(10,i)+b(11,i)
841 EEold(2,1,i)= b(12,i)-b(13,i)
842 EEold(1,2,i)= b(12,i)+b(13,i)
843 EEold(1,1,-i)= b(10,i)+b(11,i)
844 EEold(2,2,-i)=-b(10,i)+b(11,i)
845 EEold(2,1,-i)=-b(12,i)+b(13,i)
846 EEold(1,2,-i)=-b(12,i)-b(13,i)
847 write(iout,*) "TU DOCHODZE"
853 c ee(2,1,i)=ee(1,2,i)
858 &"Coefficients of the cumulants (independent of valence angles)"
859 do i=-nloctyp+1,nloctyp-1
860 write (iout,*) 'Type ',onelet(iloctyp(i))
862 write(iout,'(2f10.5)') B(3,i),B(5,i)
864 write(iout,'(2f10.5)') B(2,i),B(4,i)
867 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
871 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
875 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
880 C write (iout,*) 'KURWAKURWA'
883 C Read torsional parameters in old format
885 read (itorp,*,end=113,err=113) ntortyp,nterm_old
886 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
887 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
890 read (itorp,'(a)',end=113,err=113)
892 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
898 write (iout,'(/a/)') 'Torsional constants:'
901 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
902 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
908 C Read torsional parameters
910 IF (TOR_MODE.eq.0) THEN
912 read (itorp,*,end=113,err=113) ntortyp
913 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
915 itype2loc(i)=itortyp(i)
917 write (iout,*) 'ntortyp',ntortyp
920 itortyp(i)=-itortyp(-i)
922 c write (iout,*) 'ntortyp',ntortyp
924 do j=-ntortyp+1,ntortyp-1
925 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
927 nterm(-i,-j,iblock)=nterm(i,j,iblock)
928 nlor(-i,-j,iblock)=nlor(i,j,iblock)
931 do k=1,nterm(i,j,iblock)
932 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
934 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
935 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
936 v0ij=v0ij+si*v1(k,i,j,iblock)
939 do k=1,nlor(i,j,iblock)
940 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
941 & vlor2(k,i,j),vlor3(k,i,j)
942 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
945 v0(-i,-j,iblock)=v0ij
951 write (iout,'(/a/)') 'Torsional constants:'
955 write (iout,*) 'ityp',i,' jtyp',j," block",iblock
956 write (iout,*) 'Fourier constants'
957 do k=1,nterm(i,j,iblock)
958 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
961 write (iout,*) 'Lorenz constants'
962 do k=1,nlor(i,j,iblock)
963 write (iout,'(3(1pe15.5))')
964 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
971 C 6/23/01 Read parameters for double torsionals
975 do j=-ntortyp+1,ntortyp-1
976 do k=-ntortyp+1,ntortyp-1
977 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
978 c write (iout,*) "OK onelett",
981 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
982 & .or. t3.ne.toronelet(k)) then
983 write (iout,*) "Error in double torsional parameter file",
986 call MPI_Finalize(Ierror)
988 stop "Error in double torsional parameter file"
990 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
991 & ntermd_2(i,j,k,iblock)
992 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
993 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
994 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
995 & ntermd_1(i,j,k,iblock))
996 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
997 & ntermd_1(i,j,k,iblock))
998 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
999 & ntermd_1(i,j,k,iblock))
1000 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1001 & ntermd_1(i,j,k,iblock))
1002 C Martix of D parameters for one dimesional foureir series
1003 do l=1,ntermd_1(i,j,k,iblock)
1004 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1005 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1006 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1007 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1008 c write(iout,*) "whcodze" ,
1009 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1011 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1012 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1013 & v2s(m,l,i,j,k,iblock),
1014 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1015 C Martix of D parameters for two dimesional fourier series
1016 do l=1,ntermd_2(i,j,k,iblock)
1018 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1019 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1020 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1021 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1030 write (iout,*) 'Constants for double torsionals'
1033 do j=-ntortyp+1,ntortyp-1
1034 do k=-ntortyp+1,ntortyp-1
1035 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1036 & ' nsingle',ntermd_1(i,j,k,iblock),
1037 & ' ndouble',ntermd_2(i,j,k,iblock)
1039 write (iout,*) 'Single angles:'
1040 do l=1,ntermd_1(i,j,k,iblock)
1041 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1042 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1043 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1044 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1047 write (iout,*) 'Pairs of angles:'
1048 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1049 do l=1,ntermd_2(i,j,k,iblock)
1050 write (iout,'(i5,20f10.5)')
1051 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1054 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1055 do l=1,ntermd_2(i,j,k,iblock)
1056 write (iout,'(i5,20f10.5)')
1057 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1058 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1067 ELSE IF (TOR_MODE.eq.1) THEN
1069 C read valence-torsional parameters
1070 read (itorp,*,end=113,err=113) ntortyp
1072 write (iout,*) "Valence-torsional parameters read in ntortyp",
1074 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1075 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1077 itortyp(i)=-itortyp(-i)
1079 do i=-ntortyp+1,ntortyp-1
1080 do j=-ntortyp+1,ntortyp-1
1081 C first we read the cos and sin gamma parameters
1082 read (itorp,'(13x,a)',end=113,err=113) string
1083 write (iout,*) i,j,string
1084 read (itorp,*,end=113,err=113)
1085 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1086 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1087 do k=1,nterm_kcc(j,i)
1088 do l=1,nterm_kcc_Tb(j,i)
1089 do ll=1,nterm_kcc_Tb(j,i)
1090 read (itorp,*,end=113,err=113) ii,jj,kk,
1091 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1099 c AL 4/8/16: Calculate coefficients from one-body parameters
1102 itortyp(i)=itype2loc(i)
1105 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1107 do i=-ntortyp+1,ntortyp-1
1108 do j=-ntortyp+1,ntortyp-1
1111 do k=1,nterm_kcc_Tb(j,i)
1112 do l=1,nterm_kcc_Tb(j,i)
1113 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1114 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1115 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1116 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1119 do k=1,nterm_kcc_Tb(j,i)
1120 do l=1,nterm_kcc_Tb(j,i)
1122 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1123 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1124 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1125 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1127 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1128 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1129 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1130 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1134 cf(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
1140 if (tor_mode.gt.0 .and. lprint) then
1141 c Print valence-torsional parameters
1143 & "Parameters of the valence-torsional potentials"
1144 do i=-ntortyp+1,ntortyp-1
1145 do j=-ntortyp+1,ntortyp-1
1146 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1147 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1148 do k=1,nterm_kcc(j,i)
1149 do l=1,nterm_kcc_Tb(j,i)
1150 do ll=1,nterm_kcc_Tb(j,i)
1151 write (iout,'(3i5,2f15.4)')
1152 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1161 C Read of Side-chain backbone correlation parameters
1162 C Modified 11 May 2012 by Adasko
1165 read (isccor,*,end=119,err=119) nsccortyp
1166 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1168 isccortyp(i)=-isccortyp(-i)
1170 iscprol=isccortyp(20)
1171 c write (iout,*) 'ntortyp',ntortyp
1173 cc maxinter is maximum interaction sites
1177 read (isccor,*,end=119,err=119)
1178 &nterm_sccor(i,j),nlor_sccor(i,j)
1179 c write (iout,*) nterm_sccor(i,j)
1185 nterm_sccor(-i,j)=nterm_sccor(i,j)
1186 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1187 nterm_sccor(i,-j)=nterm_sccor(i,j)
1188 c write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1189 c & nterm_sccor(-i,-j),nterm_sccor(i,-j)
1190 do k=1,nterm_sccor(i,j)
1191 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1193 if (j.eq.iscprol) then
1194 if (i.eq.isccortyp(10)) then
1195 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1196 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1198 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1199 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1200 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1201 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1202 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1203 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1204 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1205 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1208 if (i.eq.isccortyp(10)) then
1209 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1210 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1212 if (j.eq.isccortyp(10)) then
1213 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1214 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1216 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1217 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1218 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1219 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1220 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1221 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1225 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1226 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1227 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1228 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1231 do k=1,nlor_sccor(i,j)
1232 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1233 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1234 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1235 &(1+vlor3sccor(k,i,j)**2)
1237 v0sccor(l,i,j)=v0ijsccor
1238 v0sccor(l,-i,j)=v0ijsccor1
1239 v0sccor(l,i,-j)=v0ijsccor2
1240 v0sccor(l,-i,-j)=v0ijsccor3
1246 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1250 write (iout,*) 'ityp',i,' jtyp',j
1251 write (iout,*) 'Fourier constants'
1252 do k=1,nterm_sccor(i,j)
1253 write (iout,'(2(1pe15.5))')
1254 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1256 write (iout,*) 'Lorenz constants'
1257 do k=1,nlor_sccor(i,j)
1258 write (iout,'(3(1pe15.5))')
1259 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1266 C Read electrostatic-interaction parameters
1269 write (iout,'(/a)') 'Electrostatic interaction constants:'
1270 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1271 & 'IT','JT','APP','BPP','AEL6','AEL3'
1273 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1274 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1275 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1276 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1281 app (i,j)=epp(i,j)*rri*rri
1282 bpp (i,j)=-2.0D0*epp(i,j)*rri
1283 ael6(i,j)=elpp6(i,j)*4.2D0**6
1284 ael3(i,j)=elpp3(i,j)*4.2D0**3
1285 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1286 & ael6(i,j),ael3(i,j)
1290 C Read side-chain interaction parameters.
1292 read (isidep,*,end=117,err=117) ipot,expon
1293 if (ipot.lt.1 .or. ipot.gt.5) then
1294 write (iout,'(2a)') 'Error while reading SC interaction',
1295 & 'potential file - unknown potential type.'
1299 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1300 & ', exponents are ',expon,2*expon
1301 goto (10,20,30,30,40) ipot
1302 C----------------------- LJ potential ---------------------------------
1303 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1304 & (sigma0(i),i=1,ntyp)
1306 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1307 write (iout,'(a/)') 'The epsilon array:'
1308 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1309 write (iout,'(/a)') 'One-body parameters:'
1310 write (iout,'(a,4x,a)') 'residue','sigma'
1311 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1314 C----------------------- LJK potential --------------------------------
1315 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1316 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1318 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1319 write (iout,'(a/)') 'The epsilon array:'
1320 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1321 write (iout,'(/a)') 'One-body parameters:'
1322 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1323 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1327 C---------------------- GB or BP potential -----------------------------
1329 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1331 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1332 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1333 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1334 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1336 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1337 C write(iout,*) "WARNING!!",i,ntyp
1338 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1340 C epslip(i,j)=epslip(i,j)+0.05d0
1343 C For the GB potential convert sigma'**2 into chi'
1346 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1350 write (iout,'(/a/)') 'Parameters of the BP potential:'
1351 write (iout,'(a/)') 'The epsilon array:'
1352 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1353 write (iout,'(/a)') 'One-body parameters:'
1354 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1356 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1357 & chip(i),alp(i),i=1,ntyp)
1360 C--------------------- GBV potential -----------------------------------
1361 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1362 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1363 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1365 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1366 write (iout,'(a/)') 'The epsilon array:'
1367 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1368 write (iout,'(/a)') 'One-body parameters:'
1369 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1370 & 's||/s_|_^2',' chip ',' alph '
1371 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1372 & sigii(i),chip(i),alp(i),i=1,ntyp)
1376 C-----------------------------------------------------------------------
1377 C Calculate the "working" parameters of SC interactions.
1381 epslip(i,j)=epslip(j,i)
1386 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1387 sigma(j,i)=sigma(i,j)
1388 rs0(i,j)=dwa16*sigma(i,j)
1392 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1393 & 'Working parameters of the SC interactions:',
1394 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1399 epsijlip=epslip(i,j)
1400 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1409 sigeps=dsign(1.0D0,epsij)
1411 aa_aq(i,j)=epsij*rrij*rrij
1412 bb_aq(i,j)=-sigeps*epsij*rrij
1413 aa_aq(j,i)=aa_aq(i,j)
1414 bb_aq(j,i)=bb_aq(i,j)
1415 sigeps=dsign(1.0D0,epsijlip)
1416 epsijlip=dabs(epsijlip)
1417 aa_lip(i,j)=epsijlip*rrij*rrij
1418 bb_lip(i,j)=-sigeps*epsijlip*rrij
1419 aa_lip(j,i)=aa_lip(i,j)
1420 bb_lip(j,i)=bb_lip(i,j)
1422 sigt1sq=sigma0(i)**2
1423 sigt2sq=sigma0(j)**2
1426 ratsig1=sigt2sq/sigt1sq
1427 ratsig2=1.0D0/ratsig1
1428 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1429 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1430 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1434 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1435 sigmaii(i,j)=rsum_max
1436 sigmaii(j,i)=rsum_max
1438 c sigmaii(i,j)=r0(i,j)
1439 c sigmaii(j,i)=r0(i,j)
1441 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1442 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1443 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1444 augm(i,j)=epsij*r_augm**(2*expon)
1445 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1452 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1453 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1454 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1459 C Define the SC-p interaction constants
1463 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1465 c aad(i,1)=0.3D0*4.0D0**12
1466 C Following line for constants currently implemented
1467 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1468 aad(i,1)=1.5D0*4.0D0**12
1469 c aad(i,1)=0.17D0*5.6D0**12
1471 C "Soft" SC-p repulsion
1473 C Following line for constants currently implemented
1474 c aad(i,1)=0.3D0*4.0D0**6
1475 C "Hard" SC-p repulsion
1476 bad(i,1)=3.0D0*4.0D0**6
1477 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1486 C 8/9/01 Read the SC-p interaction constants from file
1489 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1492 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1493 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1494 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1495 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1499 write (iout,*) "Parameters of SC-p interactions:"
1501 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1502 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1507 C Define the constants of the disulfide bridge
1511 c Old arbitrary potential - commented out.
1516 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1517 c energy surface of diethyl disulfide.
1518 c A. Liwo and U. Kozlowska, 11/24/03
1527 write (iout,*) dyn_ss,'dyndyn'
1529 ss_depth=ebr/wsc-0.25*eps(1,1)
1530 C write(iout,*) akcm,whpb,wsc,'KURWA'
1531 Ht=Ht/wsc-0.25*eps(1,1)
1540 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1544 write (iout,'(/a)') "Disulfide bridge parameters:"
1545 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1546 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1547 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1548 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1551 if (shield_mode.gt.0) then
1553 C VSolvSphere the volume of solving sphere
1555 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1556 C there will be no distinction between proline peptide group and normal peptide
1557 C group in case of shielding parameters
1558 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1559 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1560 write (iout,*) VSolvSphere,VSolvSphere_div
1561 C long axis of side chain
1563 long_r_sidechain(i)=vbldsc0(1,i)
1564 short_r_sidechain(i)=sigma0(i)
1569 111 write (iout,*) "Error reading bending energy parameters."
1571 112 write (iout,*) "Error reading rotamer energy parameters."
1573 113 write (iout,*) "Error reading torsional energy parameters."
1575 114 write (iout,*) "Error reading double torsional energy parameters."
1578 & "Error reading cumulant (multibody energy) parameters."
1580 116 write (iout,*) "Error reading electrostatic energy parameters."
1582 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1584 117 write (iout,*) "Error reading side chain interaction parameters."
1586 118 write (iout,*) "Error reading SCp interaction parameters."
1588 119 write (iout,*) "Error reading SCCOR parameters"
1590 121 write (iout,*) "Error reading bond parameters"
1593 call MPI_Finalize(Ierror)