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'
22 include 'COMMON.LANGEVIN'
25 character*1 onelett(4) /"G","A","P","D"/
26 character*1 toronelet(-2:2) /"p","a","G","A","P"/
28 dimension blower(3,3,maxlob)
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,mp,ip,pstok
48 read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i),
49 & msc(i),isc(i),restok(i)
54 dsc_inv(i)=1.0D0/dsc(i)
58 read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk,
61 read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
62 & aksc(j,i),abond0(j,i),j=1,nbondterm(i)),msc(i),isc(i),restok(i)
67 dsc_inv(i)=1.0D0/dsc(i)
72 write(iout,'(/a/)')"Force constants virtual bonds:"
73 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
75 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
77 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
78 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
80 write (iout,'(13x,3f10.5)')
81 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
85 read(iliptranpar,*,end=1161,err=1161) pepliptran
87 read(iliptranpar,*,end=1161,err=1161) liptranene(i)
92 C Read the parameters of the probability distribution/energy expression
93 C of the virtual-bond valence angles theta
96 read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
97 & (bthet(j,i,1,1),j=1,2)
98 read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
99 read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
100 read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
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)
115 athet(1,i,-1,-1)=athet(1,-i,1,1)
116 athet(2,i,-1,-1)=-athet(2,-i,1,1)
117 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
118 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
119 athet(1,i,-1,1)=athet(1,-i,1,1)
120 athet(2,i,-1,1)=-athet(2,-i,1,1)
121 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
122 bthet(2,i,-1,1)=bthet(2,-i,1,1)
123 athet(1,i,1,-1)=-athet(1,-i,1,1)
124 athet(2,i,1,-1)=athet(2,-i,1,1)
125 bthet(1,i,1,-1)=bthet(1,-i,1,1)
126 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
131 polthet(j,i)=polthet(j,-i)
134 gthet(j,i)=gthet(j,-i)
140 c & 'Parameters of the virtual-bond valence angles:'
141 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
142 c & ' ATHETA0 ',' A1 ',' A2 ',
145 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
146 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
148 c write (iout,'(/a/9x,5a/79(1h-))')
149 c & 'Parameters of the expression for sigma(theta_c):',
150 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
151 c & ' ALPH3 ',' SIGMA0C '
153 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
154 c & (polthet(j,i),j=0,3),sigc0(i)
156 c write (iout,'(/a/9x,5a/79(1h-))')
157 c & 'Parameters of the second gaussian:',
158 c & ' THETA0 ',' SIGMA0 ',' G1 ',
161 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
162 c & sig0(i),(gthet(j,i),j=1,3)
165 & 'Parameters of the virtual-bond valence angles:'
166 write (iout,'(/a/9x,5a/79(1h-))')
167 & 'Coefficients of expansion',
168 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
169 & ' b1*10^1 ',' b2*10^1 '
171 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
172 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
173 & (10*bthet(j,i,1,1),j=1,2)
175 write (iout,'(/a/9x,5a/79(1h-))')
176 & 'Parameters of the expression for sigma(theta_c):',
177 & ' alpha0 ',' alph1 ',' alph2 ',
178 & ' alhp3 ',' sigma0c '
180 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
181 & (polthet(j,i),j=0,3),sigc0(i)
183 write (iout,'(/a/9x,5a/79(1h-))')
184 & 'Parameters of the second gaussian:',
185 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
188 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
189 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
193 IF (tor_mode.eq.0) THEN
195 C Read the parameters of Utheta determined from ab initio surfaces
196 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
198 read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
199 & ntheterm3,nsingle,ndouble
200 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
201 read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
203 ithetyp(i)=-ithetyp(-i)
206 do i=-maxthetyp,maxthetyp
207 do j=-maxthetyp,maxthetyp
208 do k=-maxthetyp,maxthetyp
209 aa0thet(i,j,k,iblock)=0.0d0
211 aathet(l,i,j,k,iblock)=0.0d0
215 bbthet(m,l,i,j,k,iblock)=0.0d0
216 ccthet(m,l,i,j,k,iblock)=0.0d0
217 ddthet(m,l,i,j,k,iblock)=0.0d0
218 eethet(m,l,i,j,k,iblock)=0.0d0
224 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
225 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
233 C write (iout,*) "KURWA1"
236 do j=-nthetyp,nthetyp
237 do k=-nthetyp,nthetyp
238 read (ithep,'(6a)',end=111,err=111) res1
239 write(iout,*) res1,i,j,k
240 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
241 read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
243 read (ithep,*,end=111,err=111)
244 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
245 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
246 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
247 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
249 read (ithep,*,end=111,err=111)
250 & (((ffthet(llll,lll,ll,i,j,k,iblock),
251 & ffthet(lll,llll,ll,i,j,k,iblock),
252 & ggthet(llll,lll,ll,i,j,k,iblock)
253 & ,ggthet(lll,llll,ll,i,j,k,iblock),
254 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
258 C write(iout,*) "KURWA1.1"
260 C For dummy ends assign glycine-type coefficients of theta-only terms; the
261 C coefficients of theta-and-gamma-dependent terms are zero.
266 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
267 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
269 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
270 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
273 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
275 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
278 C write(iout,*) "KURWA1.5"
279 C Substitution for D aminoacids from symmetry.
282 do j=-nthetyp,nthetyp
283 do k=-nthetyp,nthetyp
284 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
286 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
290 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
291 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
292 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
293 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
299 ffthet(llll,lll,ll,i,j,k,iblock)=
300 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
301 ffthet(lll,llll,ll,i,j,k,iblock)=
302 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
303 ggthet(llll,lll,ll,i,j,k,iblock)=
304 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
305 ggthet(lll,llll,ll,i,j,k,iblock)=
306 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
316 C Control printout of the coefficients of virtual-bond-angle potentials
319 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
324 write (iout,'(//4a)')
325 & 'Type ',onelett(i),onelett(j),onelett(k)
326 write (iout,'(//a,10x,a)') " l","a[l]"
327 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
328 write (iout,'(i2,1pe15.5)')
329 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
331 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
332 & "b",l,"c",l,"d",l,"e",l
334 write (iout,'(i2,4(1pe15.5))') m,
335 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
336 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
340 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
341 & "f+",l,"f-",l,"g+",l,"g-",l
344 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
345 & ffthet(n,m,l,i,j,k,iblock),
346 & ffthet(m,n,l,i,j,k,iblock),
347 & ggthet(n,m,l,i,j,k,iblock),
348 & ggthet(m,n,l,i,j,k,iblock)
361 C here will be the apropriate recalibrating for D-aminoacid
362 read (ithep,*,end=111,err=111) nthetyp
363 do i=-nthetyp+1,nthetyp-1
364 read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
365 do j=0,nbend_kcc_Tb(i)
366 read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
371 & "Parameters of the valence-only potentials"
372 do i=-nthetyp+1,nthetyp-1
373 write (iout,'(2a)') "Type ",toronelet(i)
374 do k=0,nbend_kcc_Tb(i)
375 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
384 C write(iout,*) 'KURWA2'
387 C Read the parameters of the probability distribution/energy expression
388 C of the side chains.
391 cc write (iout,*) "tu dochodze",i
392 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
396 dsc_inv(i)=1.0D0/dsc(i)
407 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
408 & ((blower(k,l,1),l=1,k),k=1,3)
409 censc(1,1,-i)=censc(1,1,i)
410 censc(2,1,-i)=censc(2,1,i)
411 censc(3,1,-i)=-censc(3,1,i)
413 read (irotam,*,end=112,err=112) bsc(j,i)
414 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
415 & ((blower(k,l,j),l=1,k),k=1,3)
416 censc(1,j,-i)=censc(1,j,i)
417 censc(2,j,-i)=censc(2,j,i)
418 censc(3,j,-i)=-censc(3,j,i)
419 C BSC is amplitude of Gaussian
426 akl=akl+blower(k,m,j)*blower(l,m,j)
430 if (((k.eq.3).and.(l.ne.3))
431 & .or.((l.eq.3).and.(k.ne.3))) then
432 gaussc(k,l,j,-i)=-akl
433 gaussc(l,k,j,-i)=-akl
445 write (iout,'(/a)') 'Parameters of side-chain local geometry'
449 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
450 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
451 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
452 c write (iout,'(a,f10.4,4(16x,f10.4))')
453 c & 'Center ',(bsc(j,i),j=1,nlobi)
454 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
455 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
456 & 'log h',(bsc(j,i),j=1,nlobi)
457 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
458 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
465 c blower(k,l,j)=gaussc(ind,j,i)
470 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
471 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
478 C Read scrot parameters for potentials determined from all-atom AM1 calculations
479 C added by Urszula Kozlowska 07/11/2007
482 read (irotam,*,end=112,err=112)
484 read (irotam,*,end=112,err=112)
487 read(irotam,*,end=112,err=112) sc_parmin(j,i)
494 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
495 C interaction energy of the Gly, Ala, and Pro prototypes.
497 read (ifourier,*,end=115,err=115) nloctyp
498 SPLIT_FOURIERTOR = nloctyp.lt.0
499 nloctyp = iabs(nloctyp)
501 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
502 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
503 itype2loc(ntyp1)=nloctyp
504 iloctyp(nloctyp)=ntyp1
506 itype2loc(-i)=-itype2loc(i)
515 iloctyp(-i)=-iloctyp(i)
517 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
518 c write (iout,*) "nloctyp",nloctyp,
519 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
522 c write (iout,*) "NEWCORR",i
523 read (ifourier,*,end=115,err=115)
526 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
529 c write (iout,*) "NEWCORR BNEW1"
530 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
533 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
536 c write (iout,*) "NEWCORR BNEW2"
537 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
539 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
540 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
542 c write (iout,*) "NEWCORR CCNEW"
543 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
545 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
546 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
548 c write (iout,*) "NEWCORR DDNEW"
549 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
553 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
557 c write (iout,*) "NEWCORR EENEW1"
558 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
560 read (ifourier,*,end=115,err=115) e0new(ii,i)
562 c write (iout,*) (e0new(ii,i),ii=1,3)
564 c write (iout,*) "NEWCORR EENEW"
567 ccnew(ii,1,i)=ccnew(ii,1,i)/2
568 ccnew(ii,2,i)=ccnew(ii,2,i)/2
569 ddnew(ii,1,i)=ddnew(ii,1,i)/2
570 ddnew(ii,2,i)=ddnew(ii,2,i)/2
575 bnew1(ii,1,-i)= bnew1(ii,1,i)
576 bnew1(ii,2,-i)=-bnew1(ii,2,i)
577 bnew2(ii,1,-i)= bnew2(ii,1,i)
578 bnew2(ii,2,-i)=-bnew2(ii,2,i)
581 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
582 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
583 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
584 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
585 ccnew(ii,1,-i)=ccnew(ii,1,i)
586 ccnew(ii,2,-i)=-ccnew(ii,2,i)
587 ddnew(ii,1,-i)=ddnew(ii,1,i)
588 ddnew(ii,2,-i)=-ddnew(ii,2,i)
590 e0new(1,-i)= e0new(1,i)
591 e0new(2,-i)=-e0new(2,i)
592 e0new(3,-i)=-e0new(3,i)
594 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
595 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
596 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
597 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
601 write (iout,'(a)') "Coefficients of the multibody terms"
602 do i=-nloctyp+1,nloctyp-1
603 write (iout,*) "Type: ",onelet(iloctyp(i))
604 write (iout,*) "Coefficients of the expansion of B1"
606 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
608 write (iout,*) "Coefficients of the expansion of B2"
610 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
612 write (iout,*) "Coefficients of the expansion of C"
613 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
614 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
615 write (iout,*) "Coefficients of the expansion of D"
616 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
617 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
618 write (iout,*) "Coefficients of the expansion of E"
619 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
622 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
627 IF (SPLIT_FOURIERTOR) THEN
629 c write (iout,*) "NEWCORR TOR",i
630 read (ifourier,*,end=115,err=115)
633 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
636 c write (iout,*) "NEWCORR BNEW1 TOR"
637 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
640 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
643 c write (iout,*) "NEWCORR BNEW2 TOR"
644 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
646 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
647 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
649 c write (iout,*) "NEWCORR CCNEW TOR"
650 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
652 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
653 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
655 c write (iout,*) "NEWCORR DDNEW TOR"
656 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
660 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
664 c write (iout,*) "NEWCORR EENEW1 TOR"
665 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
667 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
669 c write (iout,*) (e0newtor(ii,i),ii=1,3)
671 c write (iout,*) "NEWCORR EENEW TOR"
674 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
675 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
676 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
677 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
682 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
683 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
684 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
685 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
688 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
689 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
690 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
691 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
693 e0newtor(1,-i)= e0newtor(1,i)
694 e0newtor(2,-i)=-e0newtor(2,i)
695 e0newtor(3,-i)=-e0newtor(3,i)
697 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
698 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
699 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
700 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
705 & "Single-body coefficients of the torsional potentials"
706 do i=-nloctyp+1,nloctyp-1
707 write (iout,*) "Type: ",onelet(iloctyp(i))
708 write (iout,*) "Coefficients of the expansion of B1tor"
710 write (iout,'(3hB1(,i1,1h),3f10.5)')
711 & j,(bnew1tor(k,j,i),k=1,3)
713 write (iout,*) "Coefficients of the expansion of B2tor"
715 write (iout,'(3hB2(,i1,1h),3f10.5)')
716 & j,(bnew2tor(k,j,i),k=1,3)
718 write (iout,*) "Coefficients of the expansion of Ctor"
719 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
720 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
721 write (iout,*) "Coefficients of the expansion of Dtor"
722 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
723 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
724 write (iout,*) "Coefficients of the expansion of Etor"
725 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
728 write (iout,'(1hE,2i1,2f10.5)')
729 & j,k,(eenewtor(l,j,k,i),l=1,2)
735 do i=-nloctyp+1,nloctyp-1
738 bnew1tor(ii,j,i)=bnew1(ii,j,i)
743 bnew2tor(ii,j,i)=bnew2(ii,j,i)
747 ccnewtor(ii,1,i)=ccnew(ii,1,i)
748 ccnewtor(ii,2,i)=ccnew(ii,2,i)
749 ddnewtor(ii,1,i)=ddnew(ii,1,i)
750 ddnewtor(ii,2,i)=ddnew(ii,2,i)
753 ENDIF !SPLIT_FOURIER_TOR
756 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
758 read (ifourier,*,end=115,err=115)
759 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
761 write (iout,*) 'Type ',onelet(iloctyp(i))
762 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
776 c B1tilde(1,i) = b(3)
777 c B1tilde(2,i) =-b(5)
778 c B1tilde(1,-i) =-b(3)
779 c B1tilde(2,-i) =b(5)
786 cc B1tilde(1,i) = b(3,i)
787 cc B1tilde(2,i) =-b(5,i)
788 C B1tilde(1,-i) =-b(3,i)
789 C B1tilde(2,-i) =b(5,i)
790 cc b1tilde(1,i)=0.0d0
791 cc b1tilde(2,i)=0.0d0
803 CCold(1,1,-i)= b(7,i)
804 CCold(2,2,-i)=-b(7,i)
805 CCold(2,1,-i)=-b(9,i)
806 CCold(1,2,-i)=-b(9,i)
811 c Ctilde(1,1,i)= CCold(1,1,i)
812 c Ctilde(1,2,i)= CCold(1,2,i)
813 c Ctilde(2,1,i)=-CCold(2,1,i)
814 c Ctilde(2,2,i)=-CCold(2,2,i)
816 c Ctilde(1,1,i)=0.0d0
817 c Ctilde(1,2,i)=0.0d0
818 c Ctilde(2,1,i)=0.0d0
819 c Ctilde(2,2,i)=0.0d0
824 DDold(1,1,-i)= b(6,i)
825 DDold(2,2,-i)=-b(6,i)
826 DDold(2,1,-i)=-b(8,i)
827 DDold(1,2,-i)=-b(8,i)
832 c Dtilde(1,1,i)= DD(1,1,i)
833 c Dtilde(1,2,i)= DD(1,2,i)
834 c Dtilde(2,1,i)=-DD(2,1,i)
835 c Dtilde(2,2,i)=-DD(2,2,i)
837 c Dtilde(1,1,i)=0.0d0
838 c Dtilde(1,2,i)=0.0d0
839 c Dtilde(2,1,i)=0.0d0
840 c Dtilde(2,2,i)=0.0d0
841 EEold(1,1,i)= b(10,i)+b(11,i)
842 EEold(2,2,i)=-b(10,i)+b(11,i)
843 EEold(2,1,i)= b(12,i)-b(13,i)
844 EEold(1,2,i)= b(12,i)+b(13,i)
845 EEold(1,1,-i)= b(10,i)+b(11,i)
846 EEold(2,2,-i)=-b(10,i)+b(11,i)
847 EEold(2,1,-i)=-b(12,i)+b(13,i)
848 EEold(1,2,-i)=-b(12,i)-b(13,i)
849 c write(iout,*) "TU DOCHODZE"
855 c ee(2,1,i)=ee(1,2,i)
860 &"Coefficients of the cumulants (independent of valence angles)"
861 do i=-nloctyp+1,nloctyp-1
862 write (iout,*) 'Type ',onelet(iloctyp(i))
864 write(iout,'(2f10.5)') B(3,i),B(5,i)
866 write(iout,'(2f10.5)') B(2,i),B(4,i)
869 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
873 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
877 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
882 C write (iout,*) 'KURWAKURWA'
885 C Read torsional parameters in old format
887 read (itorp,*,end=113,err=113) ntortyp,nterm_old
888 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
889 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
892 read (itorp,'(a)',end=113,err=113)
894 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
900 write (iout,'(/a/)') 'Torsional constants:'
903 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
904 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
910 C Read torsional parameters
912 IF (TOR_MODE.eq.0) THEN
914 read (itorp,*,end=113,err=113) ntortyp
915 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
917 itype2loc(i)=itortyp(i)
919 write (iout,*) 'ntortyp',ntortyp
922 itortyp(i)=-itortyp(-i)
924 c write (iout,*) 'ntortyp',ntortyp
926 do j=-ntortyp+1,ntortyp-1
927 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
929 nterm(-i,-j,iblock)=nterm(i,j,iblock)
930 nlor(-i,-j,iblock)=nlor(i,j,iblock)
933 do k=1,nterm(i,j,iblock)
934 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
936 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
937 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
938 v0ij=v0ij+si*v1(k,i,j,iblock)
941 do k=1,nlor(i,j,iblock)
942 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
943 & vlor2(k,i,j),vlor3(k,i,j)
944 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
947 v0(-i,-j,iblock)=v0ij
953 write (iout,'(/a/)') 'Torsional constants:'
957 write (iout,*) 'ityp',i,' jtyp',j," block",iblock
958 write (iout,*) 'Fourier constants'
959 do k=1,nterm(i,j,iblock)
960 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
963 write (iout,*) 'Lorenz constants'
964 do k=1,nlor(i,j,iblock)
965 write (iout,'(3(1pe15.5))')
966 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
973 C 6/23/01 Read parameters for double torsionals
977 do j=-ntortyp+1,ntortyp-1
978 do k=-ntortyp+1,ntortyp-1
979 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
980 c write (iout,*) "OK onelett",
983 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
984 & .or. t3.ne.toronelet(k)) then
985 write (iout,*) "Error in double torsional parameter file",
988 call MPI_Finalize(Ierror)
990 stop "Error in double torsional parameter file"
992 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
993 & ntermd_2(i,j,k,iblock)
994 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
995 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
996 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
997 & ntermd_1(i,j,k,iblock))
998 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
999 & ntermd_1(i,j,k,iblock))
1000 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1001 & ntermd_1(i,j,k,iblock))
1002 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1003 & ntermd_1(i,j,k,iblock))
1004 C Martix of D parameters for one dimesional foureir series
1005 do l=1,ntermd_1(i,j,k,iblock)
1006 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1007 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1008 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1009 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1010 c write(iout,*) "whcodze" ,
1011 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1013 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1014 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1015 & v2s(m,l,i,j,k,iblock),
1016 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1017 C Martix of D parameters for two dimesional fourier series
1018 do l=1,ntermd_2(i,j,k,iblock)
1020 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1021 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1022 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1023 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1032 write (iout,*) 'Constants for double torsionals'
1035 do j=-ntortyp+1,ntortyp-1
1036 do k=-ntortyp+1,ntortyp-1
1037 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1038 & ' nsingle',ntermd_1(i,j,k,iblock),
1039 & ' ndouble',ntermd_2(i,j,k,iblock)
1041 write (iout,*) 'Single angles:'
1042 do l=1,ntermd_1(i,j,k,iblock)
1043 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1044 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1045 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1046 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1049 write (iout,*) 'Pairs of angles:'
1050 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1051 do l=1,ntermd_2(i,j,k,iblock)
1052 write (iout,'(i5,20f10.5)')
1053 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1056 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1057 do l=1,ntermd_2(i,j,k,iblock)
1058 write (iout,'(i5,20f10.5)')
1059 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1060 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1069 ELSE IF (TOR_MODE.eq.1) THEN
1071 C read valence-torsional parameters
1072 read (itorp,*,end=113,err=113) ntortyp
1074 write (iout,*) "Valence-torsional parameters read in ntortyp",
1076 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1077 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1079 itortyp(i)=-itortyp(-i)
1081 do i=-ntortyp+1,ntortyp-1
1082 do j=-ntortyp+1,ntortyp-1
1083 C first we read the cos and sin gamma parameters
1084 read (itorp,'(13x,a)',end=113,err=113) string
1085 write (iout,*) i,j,string
1086 read (itorp,*,end=113,err=113)
1087 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1088 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1089 do k=1,nterm_kcc(j,i)
1090 do l=1,nterm_kcc_Tb(j,i)
1091 do ll=1,nterm_kcc_Tb(j,i)
1092 read (itorp,*,end=113,err=113) ii,jj,kk,
1093 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1101 c AL 4/8/16: Calculate coefficients from one-body parameters
1104 itortyp(i)=itype2loc(i)
1107 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1109 do i=-ntortyp+1,ntortyp-1
1110 do j=-ntortyp+1,ntortyp-1
1113 do k=1,nterm_kcc_Tb(j,i)
1114 do l=1,nterm_kcc_Tb(j,i)
1115 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1116 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1117 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1118 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1121 do k=1,nterm_kcc_Tb(j,i)
1122 do l=1,nterm_kcc_Tb(j,i)
1124 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1125 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1126 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1127 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1129 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1130 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1131 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1132 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1136 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)
1142 if (tor_mode.gt.0 .and. lprint) then
1143 c Print valence-torsional parameters
1145 & "Parameters of the valence-torsional potentials"
1146 do i=-ntortyp+1,ntortyp-1
1147 do j=-ntortyp+1,ntortyp-1
1148 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1149 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1150 do k=1,nterm_kcc(j,i)
1151 do l=1,nterm_kcc_Tb(j,i)
1152 do ll=1,nterm_kcc_Tb(j,i)
1153 write (iout,'(3i5,2f15.4)')
1154 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1163 C Read of Side-chain backbone correlation parameters
1164 C Modified 11 May 2012 by Adasko
1167 read (isccor,*,end=119,err=119) nsccortyp
1168 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1170 isccortyp(i)=-isccortyp(-i)
1172 iscprol=isccortyp(20)
1173 c write (iout,*) 'ntortyp',ntortyp
1175 cc maxinter is maximum interaction sites
1179 read (isccor,*,end=119,err=119)
1180 &nterm_sccor(i,j),nlor_sccor(i,j)
1181 c write (iout,*) nterm_sccor(i,j)
1187 nterm_sccor(-i,j)=nterm_sccor(i,j)
1188 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1189 nterm_sccor(i,-j)=nterm_sccor(i,j)
1190 c write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1191 c & nterm_sccor(-i,-j),nterm_sccor(i,-j)
1192 do k=1,nterm_sccor(i,j)
1193 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1195 if (j.eq.iscprol) then
1196 if (i.eq.isccortyp(10)) then
1197 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1198 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1200 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1201 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1202 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1203 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1204 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1205 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1206 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1207 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1210 if (i.eq.isccortyp(10)) then
1211 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1212 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1214 if (j.eq.isccortyp(10)) then
1215 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1216 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)
1222 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1223 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1227 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1228 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1229 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1230 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1233 do k=1,nlor_sccor(i,j)
1234 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1235 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1236 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1237 &(1+vlor3sccor(k,i,j)**2)
1239 v0sccor(l,i,j)=v0ijsccor
1240 v0sccor(l,-i,j)=v0ijsccor1
1241 v0sccor(l,i,-j)=v0ijsccor2
1242 v0sccor(l,-i,-j)=v0ijsccor3
1248 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1252 write (iout,*) 'ityp',i,' jtyp',j
1253 write (iout,*) 'Fourier constants'
1254 do k=1,nterm_sccor(i,j)
1255 write (iout,'(2(1pe15.5))')
1256 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1258 write (iout,*) 'Lorenz constants'
1259 do k=1,nlor_sccor(i,j)
1260 write (iout,'(3(1pe15.5))')
1261 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1268 C Read electrostatic-interaction parameters
1271 write (iout,'(/a)') 'Electrostatic interaction constants:'
1272 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1273 & 'IT','JT','APP','BPP','AEL6','AEL3'
1275 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1276 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1277 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1278 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1283 app (i,j)=epp(i,j)*rri*rri
1284 bpp (i,j)=-2.0D0*epp(i,j)*rri
1285 ael6(i,j)=elpp6(i,j)*4.2D0**6
1286 ael3(i,j)=elpp3(i,j)*4.2D0**3
1287 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1288 & ael6(i,j),ael3(i,j)
1292 C Read side-chain interaction parameters.
1294 read (isidep,*,end=117,err=117) ipot,expon
1295 if (ipot.lt.1 .or. ipot.gt.5) then
1296 write (iout,'(2a)') 'Error while reading SC interaction',
1297 & 'potential file - unknown potential type.'
1301 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1302 & ', exponents are ',expon,2*expon
1303 goto (10,20,30,30,40) ipot
1304 C----------------------- LJ potential ---------------------------------
1305 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1306 & (sigma0(i),i=1,ntyp)
1308 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1309 write (iout,'(a/)') 'The epsilon array:'
1310 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1311 write (iout,'(/a)') 'One-body parameters:'
1312 write (iout,'(a,4x,a)') 'residue','sigma'
1313 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1316 C----------------------- LJK potential --------------------------------
1317 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1318 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1320 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1321 write (iout,'(a/)') 'The epsilon array:'
1322 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1323 write (iout,'(/a)') 'One-body parameters:'
1324 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1325 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1329 C---------------------- GB or BP potential -----------------------------
1331 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1333 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1334 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1335 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1336 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1338 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1339 C write(iout,*) "WARNING!!",i,ntyp
1340 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1342 C epslip(i,j)=epslip(i,j)+0.05d0
1345 C For the GB potential convert sigma'**2 into chi'
1348 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1352 write (iout,'(/a/)') 'Parameters of the BP potential:'
1353 write (iout,'(a/)') 'The epsilon array:'
1354 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1355 write (iout,'(/a)') 'One-body parameters:'
1356 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1358 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1359 & chip(i),alp(i),i=1,ntyp)
1362 C--------------------- GBV potential -----------------------------------
1363 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1364 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1365 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1367 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1368 write (iout,'(a/)') 'The epsilon array:'
1369 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1370 write (iout,'(/a)') 'One-body parameters:'
1371 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1372 & 's||/s_|_^2',' chip ',' alph '
1373 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1374 & sigii(i),chip(i),alp(i),i=1,ntyp)
1378 C-----------------------------------------------------------------------
1379 C Calculate the "working" parameters of SC interactions.
1383 epslip(i,j)=epslip(j,i)
1388 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1389 sigma(j,i)=sigma(i,j)
1390 rs0(i,j)=dwa16*sigma(i,j)
1394 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1395 & 'Working parameters of the SC interactions:',
1396 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1401 epsijlip=epslip(i,j)
1402 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1411 sigeps=dsign(1.0D0,epsij)
1413 aa_aq(i,j)=epsij*rrij*rrij
1414 bb_aq(i,j)=-sigeps*epsij*rrij
1415 aa_aq(j,i)=aa_aq(i,j)
1416 bb_aq(j,i)=bb_aq(i,j)
1417 sigeps=dsign(1.0D0,epsijlip)
1418 epsijlip=dabs(epsijlip)
1419 aa_lip(i,j)=epsijlip*rrij*rrij
1420 bb_lip(i,j)=-sigeps*epsijlip*rrij
1421 aa_lip(j,i)=aa_lip(i,j)
1422 bb_lip(j,i)=bb_lip(i,j)
1424 sigt1sq=sigma0(i)**2
1425 sigt2sq=sigma0(j)**2
1428 ratsig1=sigt2sq/sigt1sq
1429 ratsig2=1.0D0/ratsig1
1430 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1431 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1432 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1436 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1437 sigmaii(i,j)=rsum_max
1438 sigmaii(j,i)=rsum_max
1440 c sigmaii(i,j)=r0(i,j)
1441 c sigmaii(j,i)=r0(i,j)
1443 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1444 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1445 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1446 augm(i,j)=epsij*r_augm**(2*expon)
1447 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1454 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1455 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1456 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1461 C Define the SC-p interaction constants
1465 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1467 c aad(i,1)=0.3D0*4.0D0**12
1468 C Following line for constants currently implemented
1469 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1470 aad(i,1)=1.5D0*4.0D0**12
1471 c aad(i,1)=0.17D0*5.6D0**12
1473 C "Soft" SC-p repulsion
1475 C Following line for constants currently implemented
1476 c aad(i,1)=0.3D0*4.0D0**6
1477 C "Hard" SC-p repulsion
1478 bad(i,1)=3.0D0*4.0D0**6
1479 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1488 C 8/9/01 Read the SC-p interaction constants from file
1491 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1494 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1495 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1496 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1497 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1501 write (iout,*) "Parameters of SC-p interactions:"
1503 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1504 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1509 C Define the constants of the disulfide bridge
1513 c Old arbitrary potential - commented out.
1518 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1519 c energy surface of diethyl disulfide.
1520 c A. Liwo and U. Kozlowska, 11/24/03
1529 write (iout,*) dyn_ss,'dyndyn'
1531 ss_depth=ebr/wsc-0.25*eps(1,1)
1532 C write(iout,*) akcm,whpb,wsc,'KURWA'
1533 Ht=Ht/wsc-0.25*eps(1,1)
1542 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1546 write (iout,'(/a)') "Disulfide bridge parameters:"
1547 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1548 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1549 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1550 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1553 if (shield_mode.gt.0) then
1555 C VSolvSphere the volume of solving sphere
1557 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1558 C there will be no distinction between proline peptide group and normal peptide
1559 C group in case of shielding parameters
1560 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1561 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1562 write (iout,*) VSolvSphere,VSolvSphere_div
1563 C long axis of side chain
1565 long_r_sidechain(i)=vbldsc0(1,i)
1566 short_r_sidechain(i)=sigma0(i)
1571 111 write (iout,*) "Error reading bending energy parameters."
1573 112 write (iout,*) "Error reading rotamer energy parameters."
1575 113 write (iout,*) "Error reading torsional energy parameters."
1577 114 write (iout,*) "Error reading double torsional energy parameters."
1580 & "Error reading cumulant (multibody energy) parameters."
1582 116 write (iout,*) "Error reading electrostatic energy parameters."
1584 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1586 117 write (iout,*) "Error reading side chain interaction parameters."
1588 118 write (iout,*) "Error reading SCp interaction parameters."
1590 119 write (iout,*) "Error reading SCCOR parameters"
1592 121 write (iout,*) "Error reading bond parameters"
1595 call MPI_Finalize(Ierror)