1 subroutine parmread(iparm,*)
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)
8 include 'DIMENSIONS.ZSCOPT'
9 include 'DIMENSIONS.FREE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.CHAIN'
12 include 'COMMON.INTERACT'
14 include 'COMMON.LOCAL'
15 include 'COMMON.TORSION'
16 include 'COMMON.FFIELD'
17 include 'COMMON.NAMES'
18 include 'COMMON.SBRIDGE'
19 include 'COMMON.WEIGHTS'
20 include 'COMMON.ENEPS'
21 include 'COMMON.SCCOR'
22 include 'COMMON.SCROT'
24 include 'COMMON.SHIELD'
25 include 'COMMON.CONTROL'
27 character*1 onelett(4) /"G","A","P","D"/
28 character*1 toronelet(-2:2) /"p","a","G","A","P"/
30 dimension blower(3,3,maxlob)
31 character*800 controlcard
32 character*256 bondname_t,thetname_t,rotname_t,torname_t,
33 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
39 double precision ip,mp
41 character*3 lancuch,ucase
42 C write (iout,*) "KURWA"
46 call getenv("PRINT_PARM",lancuch)
47 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
48 C Set LPRINT=.TRUE. for debugging
49 dwa16=2.0d0**(1.0d0/6.0d0)
51 C Assign virtual-bond length
55 call card_concat(controlcard,.true.)
58 key = wname(i)(:ilen(wname(i)))
59 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
62 write (iout,*) "iparm",iparm," myparm",myparm
63 c If reading not own parameters, skip assignment
64 call reada(controlcard,"D0CM",d0cm,3.78d0)
65 call reada(controlcard,"AKCM",akcm,15.1d0)
66 call reada(controlcard,"AKTH",akth,11.0d0)
67 call reada(controlcard,"AKCT",akct,12.0d0)
68 call reada(controlcard,"V1SS",v1ss,-1.08d0)
69 call reada(controlcard,"V2SS",v2ss,7.61d0)
70 call reada(controlcard,"V3SS",v3ss,13.7d0)
71 call reada(controlcard,"EBR",ebr,-5.50D0)
72 call reada(controlcard,"DTRISS",dtriss,1.0D0)
73 call reada(controlcard,"ATRISS",atriss,0.3D0)
74 call reada(controlcard,"BTRISS",btriss,0.02D0)
75 call reada(controlcard,"CTRISS",ctriss,1.0D0)
76 c dyn_ss=(index(controlcard,'DYN_SS').gt.0)
77 write (iout,*) "DYN_SS",dyn_ss
78 write(iout,*) "ATRISS",atriss
79 write(iout,*) "BTRISS",btriss
80 write(iout,*) "CTRISS",ctriss
81 write(iout,*) "DTRISS",dtriss
82 call reada(controlcard,"HT",Ht,0.0D0)
84 c Old arbitrary potential - commented out.
90 if (iparm.eq.myparm .or. .not.separate_parset) then
93 c Setup weights for UNRES
118 call card_concat(controlcard,.false.)
120 c Return if not own parameters
122 if (iparm.ne.myparm .and. separate_parset) return
124 call reads(controlcard,"BONDPAR",bondname_t,bondname)
125 open (ibond,file=bondname_t,status='old')
127 call reads(controlcard,"THETPAR",thetname_t,thetname)
128 open (ithep,file=thetname_t,status='old')
130 call reads(controlcard,"ROTPAR",rotname_t,rotname)
131 open (irotam,file=rotname_t,status='old')
133 call reads(controlcard,"TORPAR",torname_t,torname)
134 open (itorp,file=torname_t,status='old')
136 call reads(controlcard,"TORDPAR",tordname_t,tordname)
137 write (iout,*) "tor_mode",tor_mode
140 & open (itordp,file=tordname_t,status='old')
142 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
143 open (isccor,file=sccorname_t,status='old')
145 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
146 open (ifourier,file=fouriername_t,status='old')
148 call reads(controlcard,"ELEPAR",elename_t,elename)
149 open (ielep,file=elename_t,status='old')
151 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
152 open (isidep,file=sidename_t,status='old')
154 call reads(controlcard,"SCPPAR",scpname_t,scpname)
155 open (iscpp,file=scpname_t,status='old')
157 write (iout,*) "Parameter set:",iparm
158 write (iout,*) "Energy-term weights:"
160 write (iout,'(a16,f10.5)') wname(i),ww(i)
162 write (iout,*) "Sidechain potential file : ",
163 & sidename_t(:ilen(sidename_t))
165 write (iout,*) "SCp potential file : ",
166 & scpname_t(:ilen(scpname_t))
168 write (iout,*) "Electrostatic potential file : ",
169 & elename_t(:ilen(elename_t))
170 write (iout,*) "Cumulant coefficient file : ",
171 & fouriername_t(:ilen(fouriername_t))
172 write (iout,*) "Torsional parameter file : ",
173 & torname_t(:ilen(torname_t))
174 write (iout,*) "Double torsional parameter file : ",
175 & tordname_t(:ilen(tordname_t))
176 write (iout,*) "Backbone-rotamer parameter file : ",
177 & sccorname(:ilen(sccorname))
178 write (iout,*) "Bond & inertia constant file : ",
179 & bondname_t(:ilen(bondname_t))
180 write (iout,*) "Bending parameter file : ",
181 & thetname_t(:ilen(thetname_t))
182 write (iout,*) "Rotamer parameter file : ",
183 & rotname_t(:ilen(rotname_t))
186 c Read the virtual-bond parameters, masses, and moments of inertia
187 c and Stokes' radii of the peptide group and side chains
190 read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp
193 read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i)
194 dsc(i) = vbldsc0(1,i)
198 dsc_inv(i)=1.0D0/dsc(i)
202 read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk
204 read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
205 & aksc(j,i),abond0(j,i),j=1,nbondterm(i))
206 dsc(i) = vbldsc0(1,i)
210 dsc_inv(i)=1.0D0/dsc(i)
215 write(iout,'(/a/)')"Force constants virtual bonds:"
216 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
218 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
220 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
221 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
223 write (iout,'(13x,3f10.5)')
224 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
228 write (iout,*) "iliptranpar",iliptranpar
229 write (iout,*) "liptranname ",liptranname
230 read(iliptranpar,*,end=1161,err=1161) pepliptran
231 write (iout,*) "pepliptran",pepliptran
233 read(iliptranpar,*,end=1161,err=1161) liptranene(i)
234 write (iout,*) i,liptranene(i)
239 C Read the parameters of the probability distribution/energy expression
240 C of the virtual-bond valence angles theta
243 read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
244 & (bthet(j,i,1,1),j=1,2)
245 read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
246 read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
247 read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
251 athet(1,i,1,-1)=athet(1,i,1,1)
252 athet(2,i,1,-1)=athet(2,i,1,1)
253 bthet(1,i,1,-1)=-bthet(1,i,1,1)
254 bthet(2,i,1,-1)=-bthet(2,i,1,1)
255 athet(1,i,-1,1)=-athet(1,i,1,1)
256 athet(2,i,-1,1)=-athet(2,i,1,1)
257 bthet(1,i,-1,1)=bthet(1,i,1,1)
258 bthet(2,i,-1,1)=bthet(2,i,1,1)
262 athet(1,i,-1,-1)=athet(1,-i,1,1)
263 athet(2,i,-1,-1)=-athet(2,-i,1,1)
264 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
265 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
266 athet(1,i,-1,1)=athet(1,-i,1,1)
267 athet(2,i,-1,1)=-athet(2,-i,1,1)
268 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
269 bthet(2,i,-1,1)=bthet(2,-i,1,1)
270 athet(1,i,1,-1)=-athet(1,-i,1,1)
271 athet(2,i,1,-1)=athet(2,-i,1,1)
272 bthet(1,i,1,-1)=bthet(1,-i,1,1)
273 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
278 polthet(j,i)=polthet(j,-i)
281 gthet(j,i)=gthet(j,-i)
287 c & 'Parameters of the virtual-bond valence angles:'
288 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
289 c & ' ATHETA0 ',' A1 ',' A2 ',
292 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
293 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
295 c write (iout,'(/a/9x,5a/79(1h-))')
296 c & 'Parameters of the expression for sigma(theta_c):',
297 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
298 c & ' ALPH3 ',' SIGMA0C '
300 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
301 c & (polthet(j,i),j=0,3),sigc0(i)
303 c write (iout,'(/a/9x,5a/79(1h-))')
304 c & 'Parameters of the second gaussian:',
305 c & ' THETA0 ',' SIGMA0 ',' G1 ',
308 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
309 c & sig0(i),(gthet(j,i),j=1,3)
312 & 'Parameters of the virtual-bond valence angles:'
313 write (iout,'(/a/9x,5a/79(1h-))')
314 & 'Coefficients of expansion',
315 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
316 & ' b1*10^1 ',' b2*10^1 '
318 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
319 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
320 & (10*bthet(j,i,1,1),j=1,2)
322 write (iout,'(/a/9x,5a/79(1h-))')
323 & 'Parameters of the expression for sigma(theta_c):',
324 & ' alpha0 ',' alph1 ',' alph2 ',
325 & ' alhp3 ',' sigma0c '
327 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
328 & (polthet(j,i),j=0,3),sigc0(i)
330 write (iout,'(/a/9x,5a/79(1h-))')
331 & 'Parameters of the second gaussian:',
332 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
335 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
336 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
340 IF (tor_mode.eq.0) THEN
342 C Read the parameters of Utheta determined from ab initio surfaces
343 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
345 read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
346 & ntheterm3,nsingle,ndouble
347 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
348 read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
350 ithetyp(i)=-ithetyp(-i)
353 do i=-maxthetyp,maxthetyp
354 do j=-maxthetyp,maxthetyp
355 do k=-maxthetyp,maxthetyp
356 aa0thet(i,j,k,iblock)=0.0d0
358 aathet(l,i,j,k,iblock)=0.0d0
362 bbthet(m,l,i,j,k,iblock)=0.0d0
363 ccthet(m,l,i,j,k,iblock)=0.0d0
364 ddthet(m,l,i,j,k,iblock)=0.0d0
365 eethet(m,l,i,j,k,iblock)=0.0d0
371 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
372 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
380 C write (iout,*) "KURWA1"
383 do j=-nthetyp,nthetyp
384 do k=-nthetyp,nthetyp
385 read (ithep,'(6a)',end=111,err=111) res1
386 write(iout,*) res1,i,j,k
387 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
388 read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
390 read (ithep,*,end=111,err=111)
391 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
392 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
393 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
394 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
396 read (ithep,*,end=111,err=111)
397 & (((ffthet(llll,lll,ll,i,j,k,iblock),
398 & ffthet(lll,llll,ll,i,j,k,iblock),
399 & ggthet(llll,lll,ll,i,j,k,iblock)
400 & ,ggthet(lll,llll,ll,i,j,k,iblock),
401 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
405 C write(iout,*) "KURWA1.1"
407 C For dummy ends assign glycine-type coefficients of theta-only terms; the
408 C coefficients of theta-and-gamma-dependent terms are zero.
413 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
414 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
416 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
417 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
420 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
422 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
425 C write(iout,*) "KURWA1.5"
426 C Substitution for D aminoacids from symmetry.
429 do j=-nthetyp,nthetyp
430 do k=-nthetyp,nthetyp
431 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
433 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
437 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
438 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
439 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
440 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
446 ffthet(llll,lll,ll,i,j,k,iblock)=
447 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
448 ffthet(lll,llll,ll,i,j,k,iblock)=
449 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
450 ggthet(llll,lll,ll,i,j,k,iblock)=
451 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
452 ggthet(lll,llll,ll,i,j,k,iblock)=
453 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
463 C Control printout of the coefficients of virtual-bond-angle potentials
466 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
471 write (iout,'(//4a)')
472 & 'Type ',onelett(i),onelett(j),onelett(k)
473 write (iout,'(//a,10x,a)') " l","a[l]"
474 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
475 write (iout,'(i2,1pe15.5)')
476 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
478 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
479 & "b",l,"c",l,"d",l,"e",l
481 write (iout,'(i2,4(1pe15.5))') m,
482 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
483 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
487 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
488 & "f+",l,"f-",l,"g+",l,"g-",l
491 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
492 & ffthet(n,m,l,i,j,k,iblock),
493 & ffthet(m,n,l,i,j,k,iblock),
494 & ggthet(n,m,l,i,j,k,iblock),
495 & ggthet(m,n,l,i,j,k,iblock)
508 C here will be the apropriate recalibrating for D-aminoacid
509 read (ithep,*,end=111,err=111) nthetyp
510 do i=-nthetyp+1,nthetyp-1
511 read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
512 do j=0,nbend_kcc_Tb(i)
513 read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
518 & "Parameters of the valence-only potentials"
519 do i=-nthetyp+1,nthetyp-1
520 write (iout,'(2a)') "Type ",toronelet(i)
521 do k=0,nbend_kcc_Tb(i)
522 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
531 C write(iout,*) 'KURWA2'
534 C Read the parameters of the probability distribution/energy expression
535 C of the side chains.
538 cc write (iout,*) "tu dochodze",i
539 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
543 dsc_inv(i)=1.0D0/dsc(i)
554 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
555 & ((blower(k,l,1),l=1,k),k=1,3)
556 censc(1,1,-i)=censc(1,1,i)
557 censc(2,1,-i)=censc(2,1,i)
558 censc(3,1,-i)=-censc(3,1,i)
560 read (irotam,*,end=112,err=112) bsc(j,i)
561 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
562 & ((blower(k,l,j),l=1,k),k=1,3)
563 censc(1,j,-i)=censc(1,j,i)
564 censc(2,j,-i)=censc(2,j,i)
565 censc(3,j,-i)=-censc(3,j,i)
566 C BSC is amplitude of Gaussian
573 akl=akl+blower(k,m,j)*blower(l,m,j)
577 if (((k.eq.3).and.(l.ne.3))
578 & .or.((l.eq.3).and.(k.ne.3))) then
579 gaussc(k,l,j,-i)=-akl
580 gaussc(l,k,j,-i)=-akl
592 write (iout,'(/a)') 'Parameters of side-chain local geometry'
596 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
597 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
598 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
599 c write (iout,'(a,f10.4,4(16x,f10.4))')
600 c & 'Center ',(bsc(j,i),j=1,nlobi)
601 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
602 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
603 & 'log h',(bsc(j,i),j=1,nlobi)
604 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
605 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
612 c blower(k,l,j)=gaussc(ind,j,i)
617 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
618 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
625 C Read scrot parameters for potentials determined from all-atom AM1 calculations
626 C added by Urszula Kozlowska 07/11/2007
629 read (irotam,*,end=112,err=112)
631 read (irotam,*,end=112,err=112)
634 read(irotam,*,end=112,err=112) sc_parmin(j,i)
641 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
642 C interaction energy of the Gly, Ala, and Pro prototypes.
644 read (ifourier,*,end=115,err=115) nloctyp
645 SPLIT_FOURIERTOR = nloctyp.lt.0
646 nloctyp = iabs(nloctyp)
648 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
649 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
650 itype2loc(ntyp1)=nloctyp
651 iloctyp(nloctyp)=ntyp1
653 itype2loc(-i)=-itype2loc(i)
662 iloctyp(-i)=-iloctyp(i)
664 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
665 c write (iout,*) "nloctyp",nloctyp,
666 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
669 c write (iout,*) "NEWCORR",i
670 read (ifourier,*,end=115,err=115)
673 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
676 c write (iout,*) "NEWCORR BNEW1"
677 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
680 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
683 c write (iout,*) "NEWCORR BNEW2"
684 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
686 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
687 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
689 c write (iout,*) "NEWCORR CCNEW"
690 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
692 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
693 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
695 c write (iout,*) "NEWCORR DDNEW"
696 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
700 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
704 c write (iout,*) "NEWCORR EENEW1"
705 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
707 read (ifourier,*,end=115,err=115) e0new(ii,i)
709 c write (iout,*) (e0new(ii,i),ii=1,3)
711 c write (iout,*) "NEWCORR EENEW"
714 ccnew(ii,1,i)=ccnew(ii,1,i)/2
715 ccnew(ii,2,i)=ccnew(ii,2,i)/2
716 ddnew(ii,1,i)=ddnew(ii,1,i)/2
717 ddnew(ii,2,i)=ddnew(ii,2,i)/2
722 bnew1(ii,1,-i)= bnew1(ii,1,i)
723 bnew1(ii,2,-i)=-bnew1(ii,2,i)
724 bnew2(ii,1,-i)= bnew2(ii,1,i)
725 bnew2(ii,2,-i)=-bnew2(ii,2,i)
728 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
729 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
730 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
731 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
732 ccnew(ii,1,-i)=ccnew(ii,1,i)
733 ccnew(ii,2,-i)=-ccnew(ii,2,i)
734 ddnew(ii,1,-i)=ddnew(ii,1,i)
735 ddnew(ii,2,-i)=-ddnew(ii,2,i)
737 e0new(1,-i)= e0new(1,i)
738 e0new(2,-i)=-e0new(2,i)
739 e0new(3,-i)=-e0new(3,i)
741 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
742 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
743 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
744 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
748 write (iout,'(a)') "Coefficients of the multibody terms"
749 do i=-nloctyp+1,nloctyp-1
750 write (iout,*) "Type: ",onelet(iloctyp(i))
751 write (iout,*) "Coefficients of the expansion of B1"
753 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
755 write (iout,*) "Coefficients of the expansion of B2"
757 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
759 write (iout,*) "Coefficients of the expansion of C"
760 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
761 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
762 write (iout,*) "Coefficients of the expansion of D"
763 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
764 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
765 write (iout,*) "Coefficients of the expansion of E"
766 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
769 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
774 IF (SPLIT_FOURIERTOR) THEN
776 c write (iout,*) "NEWCORR TOR",i
777 read (ifourier,*,end=115,err=115)
780 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
783 c write (iout,*) "NEWCORR BNEW1 TOR"
784 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
787 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
790 c write (iout,*) "NEWCORR BNEW2 TOR"
791 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
793 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
794 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
796 c write (iout,*) "NEWCORR CCNEW TOR"
797 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
799 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
800 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
802 c write (iout,*) "NEWCORR DDNEW TOR"
803 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
807 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
811 c write (iout,*) "NEWCORR EENEW1 TOR"
812 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
814 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
816 c write (iout,*) (e0newtor(ii,i),ii=1,3)
818 c write (iout,*) "NEWCORR EENEW TOR"
821 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
822 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
823 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
824 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
829 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
830 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
831 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
832 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
835 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
836 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
837 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
838 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
840 e0newtor(1,-i)= e0newtor(1,i)
841 e0newtor(2,-i)=-e0newtor(2,i)
842 e0newtor(3,-i)=-e0newtor(3,i)
844 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
845 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
846 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
847 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
852 & "Single-body coefficients of the torsional potentials"
853 do i=-nloctyp+1,nloctyp-1
854 write (iout,*) "Type: ",onelet(iloctyp(i))
855 write (iout,*) "Coefficients of the expansion of B1tor"
857 write (iout,'(3hB1(,i1,1h),3f10.5)')
858 & j,(bnew1tor(k,j,i),k=1,3)
860 write (iout,*) "Coefficients of the expansion of B2tor"
862 write (iout,'(3hB2(,i1,1h),3f10.5)')
863 & j,(bnew2tor(k,j,i),k=1,3)
865 write (iout,*) "Coefficients of the expansion of Ctor"
866 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
867 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
868 write (iout,*) "Coefficients of the expansion of Dtor"
869 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
870 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
871 write (iout,*) "Coefficients of the expansion of Etor"
872 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
875 write (iout,'(1hE,2i1,2f10.5)')
876 & j,k,(eenewtor(l,j,k,i),l=1,2)
882 do i=-nloctyp+1,nloctyp-1
885 bnew1tor(ii,j,i)=bnew1(ii,j,i)
890 bnew2tor(ii,j,i)=bnew2(ii,j,i)
894 ccnewtor(ii,1,i)=ccnew(ii,1,i)
895 ccnewtor(ii,2,i)=ccnew(ii,2,i)
896 ddnewtor(ii,1,i)=ddnew(ii,1,i)
897 ddnewtor(ii,2,i)=ddnew(ii,2,i)
900 ENDIF !SPLIT_FOURIER_TOR
903 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
905 read (ifourier,*,end=115,err=115)
906 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
908 write (iout,*) 'Type ',onelet(iloctyp(i))
909 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
923 c B1tilde(1,i) = b(3)
924 c B1tilde(2,i) =-b(5)
925 c B1tilde(1,-i) =-b(3)
926 c B1tilde(2,-i) =b(5)
933 cc B1tilde(1,i) = b(3,i)
934 cc B1tilde(2,i) =-b(5,i)
935 C B1tilde(1,-i) =-b(3,i)
936 C B1tilde(2,-i) =b(5,i)
937 cc b1tilde(1,i)=0.0d0
938 cc b1tilde(2,i)=0.0d0
950 CCold(1,1,-i)= b(7,i)
951 CCold(2,2,-i)=-b(7,i)
952 CCold(2,1,-i)=-b(9,i)
953 CCold(1,2,-i)=-b(9,i)
958 c Ctilde(1,1,i)= CCold(1,1,i)
959 c Ctilde(1,2,i)= CCold(1,2,i)
960 c Ctilde(2,1,i)=-CCold(2,1,i)
961 c Ctilde(2,2,i)=-CCold(2,2,i)
963 c Ctilde(1,1,i)=0.0d0
964 c Ctilde(1,2,i)=0.0d0
965 c Ctilde(2,1,i)=0.0d0
966 c Ctilde(2,2,i)=0.0d0
971 DDold(1,1,-i)= b(6,i)
972 DDold(2,2,-i)=-b(6,i)
973 DDold(2,1,-i)=-b(8,i)
974 DDold(1,2,-i)=-b(8,i)
979 c Dtilde(1,1,i)= DD(1,1,i)
980 c Dtilde(1,2,i)= DD(1,2,i)
981 c Dtilde(2,1,i)=-DD(2,1,i)
982 c Dtilde(2,2,i)=-DD(2,2,i)
984 c Dtilde(1,1,i)=0.0d0
985 c Dtilde(1,2,i)=0.0d0
986 c Dtilde(2,1,i)=0.0d0
987 c Dtilde(2,2,i)=0.0d0
988 EEold(1,1,i)= b(10,i)+b(11,i)
989 EEold(2,2,i)=-b(10,i)+b(11,i)
990 EEold(2,1,i)= b(12,i)-b(13,i)
991 EEold(1,2,i)= b(12,i)+b(13,i)
992 EEold(1,1,-i)= b(10,i)+b(11,i)
993 EEold(2,2,-i)=-b(10,i)+b(11,i)
994 EEold(2,1,-i)=-b(12,i)+b(13,i)
995 EEold(1,2,-i)=-b(12,i)-b(13,i)
996 write(iout,*) "TU DOCHODZE"
1002 c ee(2,1,i)=ee(1,2,i)
1007 &"Coefficients of the cumulants (independent of valence angles)"
1008 do i=-nloctyp+1,nloctyp-1
1009 write (iout,*) 'Type ',onelet(iloctyp(i))
1011 write(iout,'(2f10.5)') B(3,i),B(5,i)
1013 write(iout,'(2f10.5)') B(2,i),B(4,i)
1016 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1020 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1024 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1029 C write (iout,*) 'KURWAKURWA'
1032 C Read torsional parameters in old format
1034 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1035 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1036 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1039 read (itorp,'(a)',end=113,err=113)
1041 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1047 write (iout,'(/a/)') 'Torsional constants:'
1050 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1051 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1057 C Read torsional parameters
1059 IF (TOR_MODE.eq.0) THEN
1061 read (itorp,*,end=113,err=113) ntortyp
1062 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1064 itype2loc(i)=itortyp(i)
1066 write (iout,*) 'ntortyp',ntortyp
1068 itype2loc(-i)=-itype2loc(i)
1070 itortyp(ntyp1)=ntortyp
1073 itortyp(i)=-itortyp(-i)
1075 c write (iout,*) 'ntortyp',ntortyp
1077 do j=-ntortyp+1,ntortyp-1
1078 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1080 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1081 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1084 do k=1,nterm(i,j,iblock)
1085 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1087 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1088 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1089 v0ij=v0ij+si*v1(k,i,j,iblock)
1092 do k=1,nlor(i,j,iblock)
1093 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1094 & vlor2(k,i,j),vlor3(k,i,j)
1095 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1098 v0(-i,-j,iblock)=v0ij
1104 write (iout,'(/a/)') 'Torsional constants:'
1107 write (iout,*) 'ityp',i,' jtyp',j
1108 write (iout,*) 'Fourier constants'
1109 do k=1,nterm(i,j,iblock)
1110 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1113 write (iout,*) 'Lorenz constants'
1114 do k=1,nlor(i,j,iblock)
1115 write (iout,'(3(1pe15.5))')
1116 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1122 C 6/23/01 Read parameters for double torsionals
1126 do j=-ntortyp+1,ntortyp-1
1127 do k=-ntortyp+1,ntortyp-1
1128 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1129 c write (iout,*) "OK onelett",
1132 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1133 & .or. t3.ne.toronelet(k)) then
1134 write (iout,*) "Error in double torsional parameter file",
1137 call MPI_Finalize(Ierror)
1139 stop "Error in double torsional parameter file"
1141 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1142 & ntermd_2(i,j,k,iblock)
1143 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1144 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1145 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1146 & ntermd_1(i,j,k,iblock))
1147 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1148 & ntermd_1(i,j,k,iblock))
1149 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1150 & ntermd_1(i,j,k,iblock))
1151 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1152 & ntermd_1(i,j,k,iblock))
1153 C Martix of D parameters for one dimesional foureir series
1154 do l=1,ntermd_1(i,j,k,iblock)
1155 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1156 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1157 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1158 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1159 c write(iout,*) "whcodze" ,
1160 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1162 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1163 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1164 & v2s(m,l,i,j,k,iblock),
1165 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1166 C Martix of D parameters for two dimesional fourier series
1167 do l=1,ntermd_2(i,j,k,iblock)
1169 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1170 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1171 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1172 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1181 write (iout,*) 'Constants for double torsionals'
1184 do j=-ntortyp+1,ntortyp-1
1185 do k=-ntortyp+1,ntortyp-1
1186 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1187 & ' nsingle',ntermd_1(i,j,k,iblock),
1188 & ' ndouble',ntermd_2(i,j,k,iblock)
1190 write (iout,*) 'Single angles:'
1191 do l=1,ntermd_1(i,j,k,iblock)
1192 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1193 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1194 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1195 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1198 write (iout,*) 'Pairs of angles:'
1199 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1200 do l=1,ntermd_2(i,j,k,iblock)
1201 write (iout,'(i5,20f10.5)')
1202 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1205 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1206 do l=1,ntermd_2(i,j,k,iblock)
1207 write (iout,'(i5,20f10.5)')
1208 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1209 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1218 ELSE IF (TOR_MODE.eq.1) THEN
1220 C read valence-torsional parameters
1221 read (itorp,*,end=113,err=113) ntortyp
1223 write (iout,*) "Valence-torsional parameters read in ntortyp",
1225 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1226 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1228 itortyp(i)=-itortyp(-i)
1230 do i=-ntortyp+1,ntortyp-1
1231 do j=-ntortyp+1,ntortyp-1
1232 C first we read the cos and sin gamma parameters
1233 read (itorp,'(13x,a)',end=113,err=113) string
1234 write (iout,*) i,j,string
1235 read (itorp,*,end=113,err=113)
1236 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1237 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1238 do k=1,nterm_kcc(j,i)
1239 do l=1,nterm_kcc_Tb(j,i)
1240 do ll=1,nterm_kcc_Tb(j,i)
1241 read (itorp,*,end=113,err=113) ii,jj,kk,
1242 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1250 c AL 4/8/16: Calculate coefficients from one-body parameters
1253 itortyp(i)=itype2loc(i)
1256 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1258 do i=-ntortyp+1,ntortyp-1
1259 do j=-ntortyp+1,ntortyp-1
1262 do k=1,nterm_kcc_Tb(j,i)
1263 do l=1,nterm_kcc_Tb(j,i)
1264 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1265 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1266 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1267 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1270 do k=1,nterm_kcc_Tb(j,i)
1271 do l=1,nterm_kcc_Tb(j,i)
1273 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1274 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1275 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1276 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1278 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1279 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1280 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1281 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1285 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)
1291 if (tor_mode.gt.0 .and. lprint) then
1292 c Print valence-torsional parameters
1294 & "Parameters of the valence-torsional potentials"
1295 do i=-ntortyp+1,ntortyp-1
1296 do j=-ntortyp+1,ntortyp-1
1297 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1298 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1299 do k=1,nterm_kcc(j,i)
1300 do l=1,nterm_kcc_Tb(j,i)
1301 do ll=1,nterm_kcc_Tb(j,i)
1302 write (iout,'(3i5,2f15.4)')
1303 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1312 C Read of Side-chain backbone correlation parameters
1313 C Modified 11 May 2012 by Adasko
1316 read (isccor,*,end=119,err=119) nsccortyp
1317 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1319 isccortyp(i)=-isccortyp(-i)
1321 iscprol=isccortyp(20)
1322 c write (iout,*) 'ntortyp',ntortyp
1324 cc maxinter is maximum interaction sites
1328 read (isccor,*,end=119,err=119)
1329 &nterm_sccor(i,j),nlor_sccor(i,j)
1330 c write (iout,*) nterm_sccor(i,j)
1336 nterm_sccor(-i,j)=nterm_sccor(i,j)
1337 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1338 nterm_sccor(i,-j)=nterm_sccor(i,j)
1339 c write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1340 c & nterm_sccor(-i,-j),nterm_sccor(i,-j)
1341 do k=1,nterm_sccor(i,j)
1342 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1344 if (j.eq.iscprol) then
1345 if (i.eq.isccortyp(10)) then
1346 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1347 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1349 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1350 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1351 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1352 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1353 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1354 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1355 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1356 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1359 if (i.eq.isccortyp(10)) then
1360 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1361 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1363 if (j.eq.isccortyp(10)) then
1364 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1365 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1367 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1368 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1369 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1370 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1371 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1372 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1376 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1377 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1378 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1379 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1382 do k=1,nlor_sccor(i,j)
1383 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1384 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1385 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1386 &(1+vlor3sccor(k,i,j)**2)
1388 v0sccor(l,i,j)=v0ijsccor
1389 v0sccor(l,-i,j)=v0ijsccor1
1390 v0sccor(l,i,-j)=v0ijsccor2
1391 v0sccor(l,-i,-j)=v0ijsccor3
1397 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1401 write (iout,*) 'ityp',i,' jtyp',j
1402 write (iout,*) 'Fourier constants'
1403 do k=1,nterm_sccor(i,j)
1404 write (iout,'(2(1pe15.5))')
1405 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1407 write (iout,*) 'Lorenz constants'
1408 do k=1,nlor_sccor(i,j)
1409 write (iout,'(3(1pe15.5))')
1410 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1417 C Read electrostatic-interaction parameters
1420 write (iout,'(/a)') 'Electrostatic interaction constants:'
1421 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1422 & 'IT','JT','APP','BPP','AEL6','AEL3'
1424 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1425 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1426 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1427 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1432 app (i,j)=epp(i,j)*rri*rri
1433 bpp (i,j)=-2.0D0*epp(i,j)*rri
1434 ael6(i,j)=elpp6(i,j)*4.2D0**6
1435 ael3(i,j)=elpp3(i,j)*4.2D0**3
1436 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1437 & ael6(i,j),ael3(i,j)
1441 C Read side-chain interaction parameters.
1443 read (isidep,*,end=117,err=117) ipot,expon
1444 if (ipot.lt.1 .or. ipot.gt.5) then
1445 write (iout,'(2a)') 'Error while reading SC interaction',
1446 & 'potential file - unknown potential type.'
1450 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1451 & ', exponents are ',expon,2*expon
1452 goto (10,20,30,30,40) ipot
1453 C----------------------- LJ potential ---------------------------------
1454 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1455 & (sigma0(i),i=1,ntyp)
1457 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1458 write (iout,'(a/)') 'The epsilon array:'
1459 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1460 write (iout,'(/a)') 'One-body parameters:'
1461 write (iout,'(a,4x,a)') 'residue','sigma'
1462 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1465 C----------------------- LJK potential --------------------------------
1466 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1467 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1469 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1470 write (iout,'(a/)') 'The epsilon array:'
1471 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1472 write (iout,'(/a)') 'One-body parameters:'
1473 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1474 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1478 C---------------------- GB or BP potential -----------------------------
1480 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1482 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1483 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1484 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1485 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1487 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1488 C write(iout,*) "WARNING!!",i,ntyp
1489 write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1491 C epslip(i,j)=epslip(i,j)+0.05d0
1494 C For the GB potential convert sigma'**2 into chi'
1497 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1501 write (iout,'(/a/)') 'Parameters of the BP potential:'
1502 write (iout,'(a/)') 'The epsilon array:'
1503 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1504 write (iout,'(/a)') 'One-body parameters:'
1505 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1507 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1508 & chip(i),alp(i),i=1,ntyp)
1511 C--------------------- GBV potential -----------------------------------
1512 c 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1513 c & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1514 c & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1516 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1518 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1519 read (isidep,*,end=117,err=117)(rr0(i),i=1,ntyp)
1520 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1521 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1522 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1524 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1525 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1528 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1529 write (iout,'(a/)') 'The epsilon array:'
1530 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1531 write (iout,'(/a)') 'One-body parameters:'
1532 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1533 & 's||/s_|_^2',' chip ',' alph '
1534 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1535 & sigii(i),chip(i),alp(i),i=1,ntyp)
1539 C-----------------------------------------------------------------------
1540 C Calculate the "working" parameters of SC interactions.
1544 epslip(i,j)=epslip(j,i)
1549 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1550 sigma(j,i)=sigma(i,j)
1551 rs0(i,j)=dwa16*sigma(i,j)
1555 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1556 & 'Working parameters of the SC interactions:',
1557 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1562 epsijlip=epslip(i,j)
1563 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1572 sigeps=dsign(1.0D0,epsij)
1574 aa_aq(i,j)=epsij*rrij*rrij
1575 bb_aq(i,j)=-sigeps*epsij*rrij
1576 aa_aq(j,i)=aa_aq(i,j)
1577 bb_aq(j,i)=bb_aq(i,j)
1578 sigeps=dsign(1.0D0,epsijlip)
1579 epsijlip=dabs(epsijlip)
1580 aa_lip(i,j)=epsijlip*rrij*rrij
1581 bb_lip(i,j)=-sigeps*epsijlip*rrij
1582 aa_lip(j,i)=aa_lip(i,j)
1583 bb_lip(j,i)=bb_lip(i,j)
1585 sigt1sq=sigma0(i)**2
1586 sigt2sq=sigma0(j)**2
1589 ratsig1=sigt2sq/sigt1sq
1590 ratsig2=1.0D0/ratsig1
1591 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1592 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1593 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1597 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1598 sigmaii(i,j)=rsum_max
1599 sigmaii(j,i)=rsum_max
1601 c sigmaii(i,j)=r0(i,j)
1602 c sigmaii(j,i)=r0(i,j)
1604 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1605 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1606 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1607 augm(i,j)=epsij*r_augm**(2*expon)
1608 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1615 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1616 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1617 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1622 C Define the SC-p interaction constants
1626 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1628 c aad(i,1)=0.3D0*4.0D0**12
1629 C Following line for constants currently implemented
1630 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1631 aad(i,1)=1.5D0*4.0D0**12
1632 c aad(i,1)=0.17D0*5.6D0**12
1634 C "Soft" SC-p repulsion
1636 C Following line for constants currently implemented
1637 c aad(i,1)=0.3D0*4.0D0**6
1638 C "Hard" SC-p repulsion
1639 bad(i,1)=3.0D0*4.0D0**6
1640 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1649 C 8/9/01 Read the SC-p interaction constants from file
1652 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1655 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1656 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1657 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1658 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1662 write (iout,*) "Parameters of SC-p interactions:"
1664 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1665 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1670 C Define the constants of the disulfide bridge
1674 c Old arbitrary potential - commented out.
1679 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1680 c energy surface of diethyl disulfide.
1681 c A. Liwo and U. Kozlowska, 11/24/03
1690 write (iout,*) dyn_ss,'Dynamic disulfide parameters'
1692 ss_depth=ebr/wsc-0.25*eps(1,1)
1693 Ht=Ht/wsc-0.25*eps(1,1)
1694 akcm=akcm*wstrain/wsc
1695 akth=akth*wstrain/wsc
1696 akct=akct*wstrain/wsc
1697 v1ss=v1ss*wstrain/wsc
1698 v2ss=v2ss*wstrain/wsc
1699 v3ss=v3ss*wstrain/wsc
1701 if (wstrain.ne.0.0) then
1702 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
1709 write (iout,'(/a)') "Disulfide bridge parameters:"
1710 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1711 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1712 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1713 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1716 if (shield_mode.gt.0) then
1718 C VSolvSphere the volume of solving sphere
1720 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1721 C there will be no distinction between proline peptide group and normal peptide
1722 C group in case of shielding parameters
1723 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1724 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1725 write (iout,*) "VSolvSphere",VSolvSphere," VSolvSphere_div",
1727 C long axis of side chain
1729 long_r_sidechain(i)=vbldsc0(1,i)
1730 short_r_sidechain(i)=sigma0(i)
1735 111 write (iout,*) "Error reading bending energy parameters."
1737 112 write (iout,*) "Error reading rotamer energy parameters."
1739 113 write (iout,*) "Error reading torsional energy parameters."
1741 114 write (iout,*) "Error reading double torsional energy parameters."
1744 & "Error reading cumulant (multibody energy) parameters."
1746 116 write (iout,*) "Error reading electrostatic energy parameters."
1748 1161 write (iout,*) "Error reading lipid parameters."
1750 117 write (iout,*) "Error reading side chain interaction parameters."
1752 118 write (iout,*) "Error reading SCp interaction parameters."
1754 119 write (iout,*) "Error reading SCCOR parameters"
1756 121 write (iout,*) "Error reading bond parameters"
1759 call MPI_Finalize(Ierror)