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 '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.WEIGHTS'
19 include 'COMMON.SCCOR'
20 include 'COMMON.SCROT'
21 include 'COMMON.SHIELD'
22 include 'COMMON.CONTROL'
23 include 'COMMON.DERIV'
24 include 'COMMON.OPTIM'
25 include 'COMMON.ENERGIES'
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
45 call getenv("PRINT_PARM",lancuch)
46 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
47 write (iout,*) "lprint ",lprint
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),0.0d0)
61 c print *,me," checkpoint 7"
63 write (iout,*) "iparm",iparm," myparm",myparm
65 c If reading not own parameters, skip assignment
66 call reada(controlcard,"D0CM",d0cm,3.78d0)
67 call reada(controlcard,"AKCM",akcm,15.1d0)
68 call reada(controlcard,"AKTH",akth,11.0d0)
69 call reada(controlcard,"AKCT",akct,12.0d0)
70 call reada(controlcard,"V1SS",v1ss,-1.08d0)
71 call reada(controlcard,"V2SS",v2ss,7.61d0)
72 call reada(controlcard,"V3SS",v3ss,13.7d0)
73 call reada(controlcard,"EBR",ebr,-5.50D0)
74 call reada(controlcard,"DTRISS",dtriss,1.0D0)
75 call reada(controlcard,"ATRISS",atriss,0.3D0)
76 call reada(controlcard,"BTRISS",btriss,0.02D0)
77 call reada(controlcard,"CTRISS",ctriss,1.0D0)
78 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
79 write(iout,*) "ATRISS",atriss
80 write(iout,*) "BTRISS",btriss
81 write(iout,*) "CTRISS",ctriss
82 write(iout,*) "DTRISS",dtriss
85 C dyn_ss_mask(i)=.false.
89 c Old arbitrary potential - commented out.
94 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
95 c energy surface of diethyl disulfide.
96 c A. Liwo and U. Kozlowska, 11/24/03
108 dyn_ssbond_ij(i,j)=1.0d300
111 call reada(controlcard,"HT",Ht,0.0D0)
113 C ss_depth=ebr/wsc-0.25*eps(1,1)
114 C write(iout,*) HT,wsc,eps(1,1),'KURWA'
115 C Ht=Ht/wsc-0.25*eps(1,1)
124 C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
127 if (iparm.eq.myparm) then
130 c Setup weights for UNRES
155 call card_concat(controlcard,.false.)
157 c Return if not own parameters
159 if (iparm.ne.myparm) return
161 call reads(controlcard,"BONDPAR",bondname_t,bondname)
162 open (ibond,file=bondname_t,status='old')
164 call reads(controlcard,"THETPAR",thetname_t,thetname)
165 open (ithep,file=thetname_t,status='old')
167 call reads(controlcard,"ROTPAR",rotname_t,rotname)
168 open (irotam,file=rotname_t,status='old')
170 call reads(controlcard,"TORPAR",torname_t,torname)
171 open (itorp,file=torname_t,status='old')
173 call reads(controlcard,"TORDPAR",tordname_t,tordname)
174 write (iout,*) "tor_mode",tor_mode
177 & open (itordp,file=tordname_t,status='old')
179 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
180 open (isccor,file=sccorname_t,status='old')
182 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
183 open (ifourier,file=fouriername_t,status='old')
185 call reads(controlcard,"ELEPAR",elename_t,elename)
186 open (ielep,file=elename_t,status='old')
188 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
189 open (isidep,file=sidename_t,status='old')
191 call reads(controlcard,"SCPPAR",scpname_t,scpname)
192 open (iscpp,file=scpname_t,status='old')
194 write (iout,*) "Parameter set:",iparm
195 write (iout,*) "Energy-term weights:"
197 write (iout,'(a16,f10.5)') wname(i),ww(i)
199 write (iout,*) "Sidechain potential file : ",
200 & sidename_t(:ilen(sidename_t))
202 write (iout,*) "SCp potential file : ",
203 & scpname_t(:ilen(scpname_t))
205 write (iout,*) "Electrostatic potential file : ",
206 & elename_t(:ilen(elename_t))
207 write (iout,*) "Cumulant coefficient file : ",
208 & fouriername_t(:ilen(fouriername_t))
209 write (iout,*) "Torsional parameter file : ",
210 & torname_t(:ilen(torname_t))
211 write (iout,*) "Double torsional parameter file : ",
212 & tordname_t(:ilen(tordname_t))
213 write (iout,*) "Backbone-rotamer parameter file : ",
214 & sccorname(:ilen(sccorname))
215 write (iout,*) "Bond & inertia constant file : ",
216 & bondname_t(:ilen(bondname_t))
217 write (iout,*) "Bending parameter file : ",
218 & thetname_t(:ilen(thetname_t))
219 write (iout,*) "Rotamer parameter file : ",
220 & rotname_t(:ilen(rotname_t))
223 c Read the virtual-bond parameters, masses, and moments of inertia
224 c and Stokes' radii of the peptide group and side chains
227 read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp
230 read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i)
231 dsc(i) = vbldsc0(1,i)
235 dsc_inv(i)=1.0D0/dsc(i)
239 read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk
241 read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
242 & aksc(j,i),abond0(j,i),j=1,nbondterm(i))
243 dsc(i) = vbldsc0(1,i)
247 dsc_inv(i)=1.0D0/dsc(i)
252 write(iout,'(/a/)')"Force constants virtual bonds:"
253 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
255 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
257 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
258 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
260 write (iout,'(13x,3f10.5)')
261 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
265 read(iliptranpar,*,end=1161,err=1161) pepliptran
267 read(iliptranpar,*,end=1161,err=1161) liptranene(i)
272 C Read the parameters of the probability distribution/energy expression
273 C of the virtual-bond valence angles theta
276 read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
277 & (bthet(j,i,1,1),j=1,2)
278 read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
279 read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
280 read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
284 athet(1,i,1,-1)=athet(1,i,1,1)
285 athet(2,i,1,-1)=athet(2,i,1,1)
286 bthet(1,i,1,-1)=-bthet(1,i,1,1)
287 bthet(2,i,1,-1)=-bthet(2,i,1,1)
288 athet(1,i,-1,1)=-athet(1,i,1,1)
289 athet(2,i,-1,1)=-athet(2,i,1,1)
290 bthet(1,i,-1,1)=bthet(1,i,1,1)
291 bthet(2,i,-1,1)=bthet(2,i,1,1)
295 athet(1,i,-1,-1)=athet(1,-i,1,1)
296 athet(2,i,-1,-1)=-athet(2,-i,1,1)
297 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
298 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
299 athet(1,i,-1,1)=athet(1,-i,1,1)
300 athet(2,i,-1,1)=-athet(2,-i,1,1)
301 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
302 bthet(2,i,-1,1)=bthet(2,-i,1,1)
303 athet(1,i,1,-1)=-athet(1,-i,1,1)
304 athet(2,i,1,-1)=athet(2,-i,1,1)
305 bthet(1,i,1,-1)=bthet(1,-i,1,1)
306 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
311 polthet(j,i)=polthet(j,-i)
314 gthet(j,i)=gthet(j,-i)
320 c & 'Parameters of the virtual-bond valence angles:'
321 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
322 c & ' ATHETA0 ',' A1 ',' A2 ',
325 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
326 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
328 c write (iout,'(/a/9x,5a/79(1h-))')
329 c & 'Parameters of the expression for sigma(theta_c):',
330 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
331 c & ' ALPH3 ',' SIGMA0C '
333 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
334 c & (polthet(j,i),j=0,3),sigc0(i)
336 c write (iout,'(/a/9x,5a/79(1h-))')
337 c & 'Parameters of the second gaussian:',
338 c & ' THETA0 ',' SIGMA0 ',' G1 ',
341 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
342 c & sig0(i),(gthet(j,i),j=1,3)
345 & 'Parameters of the virtual-bond valence angles:'
346 write (iout,'(/a/9x,5a/79(1h-))')
347 & 'Coefficients of expansion',
348 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
349 & ' b1*10^1 ',' b2*10^1 '
351 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
352 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
353 & (10*bthet(j,i,1,1),j=1,2)
355 write (iout,'(/a/9x,5a/79(1h-))')
356 & 'Parameters of the expression for sigma(theta_c):',
357 & ' alpha0 ',' alph1 ',' alph2 ',
358 & ' alhp3 ',' sigma0c '
360 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
361 & (polthet(j,i),j=0,3),sigc0(i)
363 write (iout,'(/a/9x,5a/79(1h-))')
364 & 'Parameters of the second gaussian:',
365 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
368 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
369 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
373 IF (tor_mode.eq.0) THEN
375 C Read the parameters of Utheta determined from ab initio surfaces
376 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
378 read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
379 & ntheterm3,nsingle,ndouble
380 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
381 read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
383 ithetyp(i)=-ithetyp(-i)
386 do i=-maxthetyp,maxthetyp
387 do j=-maxthetyp,maxthetyp
388 do k=-maxthetyp,maxthetyp
389 aa0thet(i,j,k,iblock)=0.0d0
391 aathet(l,i,j,k,iblock)=0.0d0
395 bbthet(m,l,i,j,k,iblock)=0.0d0
396 ccthet(m,l,i,j,k,iblock)=0.0d0
397 ddthet(m,l,i,j,k,iblock)=0.0d0
398 eethet(m,l,i,j,k,iblock)=0.0d0
404 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
405 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
413 C write (iout,*) "KURWA1"
416 do j=-nthetyp,nthetyp
417 do k=-nthetyp,nthetyp
418 read (ithep,'(6a)',end=111,err=111) res1
419 write(iout,*) res1,i,j,k
420 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
421 read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
423 read (ithep,*,end=111,err=111)
424 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
425 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
426 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
427 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
429 read (ithep,*,end=111,err=111)
430 & (((ffthet(llll,lll,ll,i,j,k,iblock),
431 & ffthet(lll,llll,ll,i,j,k,iblock),
432 & ggthet(llll,lll,ll,i,j,k,iblock)
433 & ,ggthet(lll,llll,ll,i,j,k,iblock),
434 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
438 C write(iout,*) "KURWA1.1"
440 C For dummy ends assign glycine-type coefficients of theta-only terms; the
441 C coefficients of theta-and-gamma-dependent terms are zero.
446 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
447 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
449 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
450 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
453 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
455 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
458 C write(iout,*) "KURWA1.5"
459 C Substitution for D aminoacids from symmetry.
462 do j=-nthetyp,nthetyp
463 do k=-nthetyp,nthetyp
464 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
466 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
470 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
471 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
472 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
473 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
479 ffthet(llll,lll,ll,i,j,k,iblock)=
480 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
481 ffthet(lll,llll,ll,i,j,k,iblock)=
482 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
483 ggthet(llll,lll,ll,i,j,k,iblock)=
484 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
485 ggthet(lll,llll,ll,i,j,k,iblock)=
486 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
496 C Control printout of the coefficients of virtual-bond-angle potentials
499 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
504 write (iout,'(//4a)')
505 & 'Type ',onelett(i),onelett(j),onelett(k)
506 write (iout,'(//a,10x,a)') " l","a[l]"
507 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
508 write (iout,'(i2,1pe15.5)')
509 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
511 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
512 & "b",l,"c",l,"d",l,"e",l
514 write (iout,'(i2,4(1pe15.5))') m,
515 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
516 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
520 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
521 & "f+",l,"f-",l,"g+",l,"g-",l
524 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
525 & ffthet(n,m,l,i,j,k,iblock),
526 & ffthet(m,n,l,i,j,k,iblock),
527 & ggthet(n,m,l,i,j,k,iblock),
528 & ggthet(m,n,l,i,j,k,iblock)
541 C here will be the apropriate recalibrating for D-aminoacid
542 read (ithep,*,end=111,err=111) nthetyp
543 do i=-nthetyp+1,nthetyp-1
544 read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
545 do j=0,nbend_kcc_Tb(i)
546 read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
551 & "Parameters of the valence-only potentials"
552 do i=-nthetyp+1,nthetyp-1
553 write (iout,'(2a)') "Type ",toronelet(i)
554 do k=0,nbend_kcc_Tb(i)
555 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
564 C write(iout,*) 'KURWA2'
567 C Read the parameters of the probability distribution/energy expression
568 C of the side chains.
571 cc write (iout,*) "tu dochodze",i
572 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
576 dsc_inv(i)=1.0D0/dsc(i)
587 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
588 & ((blower(k,l,1),l=1,k),k=1,3)
589 censc(1,1,-i)=censc(1,1,i)
590 censc(2,1,-i)=censc(2,1,i)
591 censc(3,1,-i)=-censc(3,1,i)
593 read (irotam,*,end=112,err=112) bsc(j,i)
594 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
595 & ((blower(k,l,j),l=1,k),k=1,3)
596 censc(1,j,-i)=censc(1,j,i)
597 censc(2,j,-i)=censc(2,j,i)
598 censc(3,j,-i)=-censc(3,j,i)
599 C BSC is amplitude of Gaussian
606 akl=akl+blower(k,m,j)*blower(l,m,j)
610 if (((k.eq.3).and.(l.ne.3))
611 & .or.((l.eq.3).and.(k.ne.3))) then
612 gaussc(k,l,j,-i)=-akl
613 gaussc(l,k,j,-i)=-akl
625 write (iout,'(/a)') 'Parameters of side-chain local geometry'
629 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
630 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
631 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
632 c write (iout,'(a,f10.4,4(16x,f10.4))')
633 c & 'Center ',(bsc(j,i),j=1,nlobi)
634 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
635 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
636 & 'log h',(bsc(j,i),j=1,nlobi)
637 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
638 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
645 c blower(k,l,j)=gaussc(ind,j,i)
650 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
651 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
658 C Read scrot parameters for potentials determined from all-atom AM1 calculations
659 C added by Urszula Kozlowska 07/11/2007
662 read (irotam,*,end=112,err=112)
664 read (irotam,*,end=112,err=112)
667 read(irotam,*,end=112,err=112) sc_parmin(j,i)
674 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
675 C interaction energy of the Gly, Ala, and Pro prototypes.
677 read (ifourier,*,end=115,err=115) nloctyp
678 SPLIT_FOURIERTOR = nloctyp.lt.0
679 nloctyp = iabs(nloctyp)
680 write (iout,*) "nloctyp",nloctyp
682 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
683 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
684 itype2loc(ntyp1)=nloctyp
685 iloctyp(nloctyp)=ntyp1
687 itype2loc(-i)=-itype2loc(i)
696 iloctyp(-i)=-iloctyp(i)
698 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
699 c write (iout,*) "nloctyp",nloctyp,
700 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
703 c write (iout,*) "NEWCORR",i
704 read (ifourier,*,end=115,err=115)
707 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
710 c write (iout,*) "NEWCORR BNEW1"
711 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
714 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
717 c write (iout,*) "NEWCORR BNEW2"
718 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
720 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
721 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
723 c write (iout,*) "NEWCORR CCNEW"
724 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
726 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
727 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
729 c write (iout,*) "NEWCORR DDNEW"
730 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
734 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
738 c write (iout,*) "NEWCORR EENEW1"
739 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
741 read (ifourier,*,end=115,err=115) e0new(ii,i)
743 c write (iout,*) (e0new(ii,i),ii=1,3)
745 c write (iout,*) "NEWCORR EENEW"
748 ccnew(ii,1,i)=ccnew(ii,1,i)/2
749 ccnew(ii,2,i)=ccnew(ii,2,i)/2
750 ddnew(ii,1,i)=ddnew(ii,1,i)/2
751 ddnew(ii,2,i)=ddnew(ii,2,i)/2
756 bnew1(ii,1,-i)= bnew1(ii,1,i)
757 bnew1(ii,2,-i)=-bnew1(ii,2,i)
758 bnew2(ii,1,-i)= bnew2(ii,1,i)
759 bnew2(ii,2,-i)=-bnew2(ii,2,i)
762 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
763 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
764 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
765 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
766 ccnew(ii,1,-i)=ccnew(ii,1,i)
767 ccnew(ii,2,-i)=-ccnew(ii,2,i)
768 ddnew(ii,1,-i)=ddnew(ii,1,i)
769 ddnew(ii,2,-i)=-ddnew(ii,2,i)
771 e0new(1,-i)= e0new(1,i)
772 e0new(2,-i)=-e0new(2,i)
773 e0new(3,-i)=-e0new(3,i)
775 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
776 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
777 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
778 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
782 write (iout,'(a)') "Coefficients of the multibody terms"
783 do i=-nloctyp+1,nloctyp-1
784 write (iout,*) "Type: ",onelet(iloctyp(i))
785 write (iout,*) "Coefficients of the expansion of B1"
787 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
789 write (iout,*) "Coefficients of the expansion of B2"
791 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
793 write (iout,*) "Coefficients of the expansion of C"
794 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
795 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
796 write (iout,*) "Coefficients of the expansion of D"
797 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
798 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
799 write (iout,*) "Coefficients of the expansion of E"
800 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
803 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
808 IF (SPLIT_FOURIERTOR) THEN
810 c write (iout,*) "NEWCORR TOR",i
811 read (ifourier,*,end=115,err=115)
814 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
817 c write (iout,*) "NEWCORR BNEW1 TOR"
818 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
821 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
824 c write (iout,*) "NEWCORR BNEW2 TOR"
825 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
827 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
828 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
830 c write (iout,*) "NEWCORR CCNEW TOR"
831 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
833 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
834 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
836 c write (iout,*) "NEWCORR DDNEW TOR"
837 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
841 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
845 c write (iout,*) "NEWCORR EENEW1 TOR"
846 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
848 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
850 c write (iout,*) (e0newtor(ii,i),ii=1,3)
852 c write (iout,*) "NEWCORR EENEW TOR"
855 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
856 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
857 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
858 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
863 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
864 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
865 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
866 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
869 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
870 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
871 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
872 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
874 e0newtor(1,-i)= e0newtor(1,i)
875 e0newtor(2,-i)=-e0newtor(2,i)
876 e0newtor(3,-i)=-e0newtor(3,i)
878 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
879 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
880 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
881 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
886 & "Single-body coefficients of the torsional potentials"
887 do i=-nloctyp+1,nloctyp-1
888 write (iout,*) "Type: ",onelet(iloctyp(i))
889 write (iout,*) "Coefficients of the expansion of B1tor"
891 write (iout,'(3hB1(,i1,1h),3f10.5)')
892 & j,(bnew1tor(k,j,i),k=1,3)
894 write (iout,*) "Coefficients of the expansion of B2tor"
896 write (iout,'(3hB2(,i1,1h),3f10.5)')
897 & j,(bnew2tor(k,j,i),k=1,3)
899 write (iout,*) "Coefficients of the expansion of Ctor"
900 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
901 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
902 write (iout,*) "Coefficients of the expansion of Dtor"
903 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
904 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
905 write (iout,*) "Coefficients of the expansion of Etor"
906 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
909 write (iout,'(1hE,2i1,2f10.5)')
910 & j,k,(eenewtor(l,j,k,i),l=1,2)
916 do i=-nloctyp+1,nloctyp-1
919 bnew1tor(ii,j,i)=bnew1(ii,j,i)
924 bnew2tor(ii,j,i)=bnew2(ii,j,i)
928 ccnewtor(ii,1,i)=ccnew(ii,1,i)
929 ccnewtor(ii,2,i)=ccnew(ii,2,i)
930 ddnewtor(ii,1,i)=ddnew(ii,1,i)
931 ddnewtor(ii,2,i)=ddnew(ii,2,i)
934 ENDIF !SPLIT_FOURIER_TOR
937 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
939 read (ifourier,*,end=115,err=115)
940 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
942 write (iout,*) 'Type ',onelet(iloctyp(i))
943 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
957 c B1tilde(1,i) = b(3)
958 c B1tilde(2,i) =-b(5)
959 c B1tilde(1,-i) =-b(3)
960 c B1tilde(2,-i) =b(5)
967 cc B1tilde(1,i) = b(3,i)
968 cc B1tilde(2,i) =-b(5,i)
969 C B1tilde(1,-i) =-b(3,i)
970 C B1tilde(2,-i) =b(5,i)
971 cc b1tilde(1,i)=0.0d0
972 cc b1tilde(2,i)=0.0d0
984 CCold(1,1,-i)= b(7,i)
985 CCold(2,2,-i)=-b(7,i)
986 CCold(2,1,-i)=-b(9,i)
987 CCold(1,2,-i)=-b(9,i)
992 c Ctilde(1,1,i)= CCold(1,1,i)
993 c Ctilde(1,2,i)= CCold(1,2,i)
994 c Ctilde(2,1,i)=-CCold(2,1,i)
995 c Ctilde(2,2,i)=-CCold(2,2,i)
997 c Ctilde(1,1,i)=0.0d0
998 c Ctilde(1,2,i)=0.0d0
999 c Ctilde(2,1,i)=0.0d0
1000 c Ctilde(2,2,i)=0.0d0
1001 DDold(1,1,i)= b(6,i)
1002 DDold(2,2,i)=-b(6,i)
1003 DDold(2,1,i)= b(8,i)
1004 DDold(1,2,i)= b(8,i)
1005 DDold(1,1,-i)= b(6,i)
1006 DDold(2,2,-i)=-b(6,i)
1007 DDold(2,1,-i)=-b(8,i)
1008 DDold(1,2,-i)=-b(8,i)
1013 c Dtilde(1,1,i)= DD(1,1,i)
1014 c Dtilde(1,2,i)= DD(1,2,i)
1015 c Dtilde(2,1,i)=-DD(2,1,i)
1016 c Dtilde(2,2,i)=-DD(2,2,i)
1018 c Dtilde(1,1,i)=0.0d0
1019 c Dtilde(1,2,i)=0.0d0
1020 c Dtilde(2,1,i)=0.0d0
1021 c Dtilde(2,2,i)=0.0d0
1022 EEold(1,1,i)= b(10,i)+b(11,i)
1023 EEold(2,2,i)=-b(10,i)+b(11,i)
1024 EEold(2,1,i)= b(12,i)-b(13,i)
1025 EEold(1,2,i)= b(12,i)+b(13,i)
1026 EEold(1,1,-i)= b(10,i)+b(11,i)
1027 EEold(2,2,-i)=-b(10,i)+b(11,i)
1028 EEold(2,1,-i)=-b(12,i)+b(13,i)
1029 EEold(1,2,-i)=-b(12,i)-b(13,i)
1030 write(iout,*) "TU DOCHODZE"
1036 c ee(2,1,i)=ee(1,2,i)
1041 &"Coefficients of the cumulants (independent of valence angles)"
1042 do i=-nloctyp+1,nloctyp-1
1043 write (iout,*) 'Type ',onelet(iloctyp(i))
1045 write(iout,'(2f10.5)') B(3,i),B(5,i)
1047 write(iout,'(2f10.5)') B(2,i),B(4,i)
1050 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1054 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1058 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1063 C write (iout,*) 'KURWAKURWA'
1066 C Read torsional parameters in old format
1068 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1069 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1070 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1073 read (itorp,'(a)',end=113,err=113)
1075 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1081 write (iout,'(/a/)') 'Torsional constants:'
1084 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1085 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1091 C Read torsional parameters
1093 IF (TOR_MODE.eq.0) THEN
1095 read (itorp,*,end=113,err=113) ntortyp
1096 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1097 write (iout,*) 'ntortyp',ntortyp
1099 itype2loc(i)=itortyp(i)
1102 itype2loc(-i)=-itype2loc(i)
1106 itortyp(i)=-itortyp(-i)
1108 c write (iout,*) 'ntortyp',ntortyp
1110 do j=-ntortyp+1,ntortyp-1
1111 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1113 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1114 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1117 do k=1,nterm(i,j,iblock)
1118 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1120 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1121 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1122 v0ij=v0ij+si*v1(k,i,j,iblock)
1125 do k=1,nlor(i,j,iblock)
1126 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1127 & vlor2(k,i,j),vlor3(k,i,j)
1128 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1131 v0(-i,-j,iblock)=v0ij
1137 write (iout,'(/a/)') 'Torsional constants:'
1141 write (iout,*) 'ityp',i,' jtyp',j
1142 write (iout,*) 'Fourier constants'
1143 do k=1,nterm(i,j,iblock)
1144 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1147 write (iout,*) 'Lorenz constants'
1148 do k=1,nlor(i,j,iblock)
1149 write (iout,'(3(1pe15.5))')
1150 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1157 C 6/23/01 Read parameters for double torsionals
1161 do j=-ntortyp+1,ntortyp-1
1162 do k=-ntortyp+1,ntortyp-1
1163 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1164 c write (iout,*) "OK onelett",
1167 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1168 & .or. t3.ne.toronelet(k)) then
1169 write (iout,*) "Error in double torsional parameter file",
1172 call MPI_Finalize(Ierror)
1174 stop "Error in double torsional parameter file"
1176 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1177 & ntermd_2(i,j,k,iblock)
1178 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1179 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1180 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1181 & ntermd_1(i,j,k,iblock))
1182 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1183 & ntermd_1(i,j,k,iblock))
1184 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1185 & ntermd_1(i,j,k,iblock))
1186 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1187 & ntermd_1(i,j,k,iblock))
1188 C Martix of D parameters for one dimesional foureir series
1189 do l=1,ntermd_1(i,j,k,iblock)
1190 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1191 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1192 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1193 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1194 c write(iout,*) "whcodze" ,
1195 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1197 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1198 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1199 & v2s(m,l,i,j,k,iblock),
1200 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1201 C Martix of D parameters for two dimesional fourier series
1202 do l=1,ntermd_2(i,j,k,iblock)
1204 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1205 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1206 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1207 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1216 write (iout,*) 'Constants for double torsionals'
1219 do j=-ntortyp+1,ntortyp-1
1220 do k=-ntortyp+1,ntortyp-1
1221 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1222 & ' nsingle',ntermd_1(i,j,k,iblock),
1223 & ' ndouble',ntermd_2(i,j,k,iblock)
1225 write (iout,*) 'Single angles:'
1226 do l=1,ntermd_1(i,j,k,iblock)
1227 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1228 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1229 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1230 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1233 write (iout,*) 'Pairs of angles:'
1234 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1235 do l=1,ntermd_2(i,j,k,iblock)
1236 write (iout,'(i5,20f10.5)')
1237 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1240 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1241 do l=1,ntermd_2(i,j,k,iblock)
1242 write (iout,'(i5,20f10.5)')
1243 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1244 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1253 ELSE IF (TOR_MODE.eq.1) THEN
1255 C read valence-torsional parameters
1256 read (itorp,*,end=113,err=113) ntortyp
1258 write (iout,*) "Valence-torsional parameters read in ntortyp",
1260 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1261 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1263 itortyp(i)=-itortyp(-i)
1265 do i=-ntortyp+1,ntortyp-1
1266 do j=-ntortyp+1,ntortyp-1
1267 C first we read the cos and sin gamma parameters
1268 read (itorp,'(13x,a)',end=113,err=113) string
1269 write (iout,*) i,j,string
1270 read (itorp,*,end=113,err=113)
1271 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1272 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1273 do k=1,nterm_kcc(j,i)
1274 do l=1,nterm_kcc_Tb(j,i)
1275 do ll=1,nterm_kcc_Tb(j,i)
1276 read (itorp,*,end=113,err=113) ii,jj,kk,
1277 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1285 c AL 4/8/16: Calculate coefficients from one-body parameters
1288 itortyp(i)=itype2loc(i)
1292 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1294 do i=-ntortyp+1,ntortyp-1
1295 do j=-ntortyp+1,ntortyp-1
1298 do k=1,nterm_kcc_Tb(j,i)
1299 do l=1,nterm_kcc_Tb(j,i)
1300 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1301 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1302 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1303 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1306 do k=1,nterm_kcc_Tb(j,i)
1307 do l=1,nterm_kcc_Tb(j,i)
1308 c v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1309 c & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1310 c v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1311 c & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1312 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1313 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1314 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1315 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1318 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)
1324 if (tor_mode.gt.0 .and. lprint) then
1325 c Print valence-torsional parameters
1327 & "Parameters of the valence-torsional potentials"
1328 do i=-ntortyp+1,ntortyp-1
1329 do j=-ntortyp+1,ntortyp-1
1330 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1331 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1332 do k=1,nterm_kcc(j,i)
1333 do l=1,nterm_kcc_Tb(j,i)
1334 do ll=1,nterm_kcc_Tb(j,i)
1335 write (iout,'(3i5,2f15.4)')
1336 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1345 C Read of Side-chain backbone correlation parameters
1346 C Modified 11 May 2012 by Adasko
1349 read (isccor,*,end=119,err=119) nsccortyp
1350 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1352 isccortyp(i)=-isccortyp(-i)
1354 iscprol=isccortyp(20)
1355 c write (iout,*) 'ntortyp',ntortyp
1357 cc maxinter is maximum interaction sites
1361 read (isccor,*,end=119,err=119)
1362 &nterm_sccor(i,j),nlor_sccor(i,j)
1363 c write (iout,*) nterm_sccor(i,j)
1369 nterm_sccor(-i,j)=nterm_sccor(i,j)
1370 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1371 nterm_sccor(i,-j)=nterm_sccor(i,j)
1372 c write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1373 c & nterm_sccor(-i,-j),nterm_sccor(i,-j)
1374 do k=1,nterm_sccor(i,j)
1375 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1377 if (j.eq.iscprol) then
1378 if (i.eq.isccortyp(10)) then
1379 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1380 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1382 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1383 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1384 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1385 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1386 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1387 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1388 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1389 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1392 if (i.eq.isccortyp(10)) then
1393 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1394 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1396 if (j.eq.isccortyp(10)) then
1397 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1398 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1400 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1401 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1402 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1403 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1404 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1405 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1409 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1410 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1411 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1412 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1415 do k=1,nlor_sccor(i,j)
1416 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1417 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1418 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1419 &(1+vlor3sccor(k,i,j)**2)
1421 v0sccor(l,i,j)=v0ijsccor
1422 v0sccor(l,-i,j)=v0ijsccor1
1423 v0sccor(l,i,-j)=v0ijsccor2
1424 v0sccor(l,-i,-j)=v0ijsccor3
1430 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1434 write (iout,*) 'ityp',i,' jtyp',j
1435 write (iout,*) 'Fourier constants'
1436 do k=1,nterm_sccor(i,j)
1437 write (iout,'(2(1pe15.5))')
1438 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1440 write (iout,*) 'Lorenz constants'
1441 do k=1,nlor_sccor(i,j)
1442 write (iout,'(3(1pe15.5))')
1443 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1450 C Read electrostatic-interaction parameters
1453 write (iout,'(/a)') 'Electrostatic interaction constants:'
1454 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1455 & 'IT','JT','APP','BPP','AEL6','AEL3'
1457 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1458 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1459 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1460 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1465 app (i,j)=epp(i,j)*rri*rri
1466 bpp (i,j)=-2.0D0*epp(i,j)*rri
1467 ael6(i,j)=elpp6(i,j)*4.2D0**6
1468 ael3(i,j)=elpp3(i,j)*4.2D0**3
1470 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1471 & ael6(i,j),ael3(i,j)
1475 C Read side-chain interaction parameters.
1477 read (isidep,*) ipot,expon
1478 if (ipot.lt.1 .or. ipot.gt.6) then
1479 write (iout,'(2a)') 'Error while reading SC interaction',
1480 & 'potential file - unknown potential type.'
1484 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1485 & ', exponents are ',expon,2*expon
1486 goto (10,20,30,30,40,50,60) ipot
1487 C----------------------- LJ potential ---------------------------------
1488 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1490 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1491 write (iout,'(a/)') 'The epsilon array:'
1492 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1493 write (iout,'(/a)') 'One-body parameters:'
1494 write (iout,'(a,4x,a)') 'residue','sigma'
1495 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1498 C----------------------- LJK potential --------------------------------
1499 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1500 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1502 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1503 write (iout,'(a/)') 'The epsilon array:'
1504 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1505 write (iout,'(/a)') 'One-body parameters:'
1506 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1507 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1511 C---------------------- GB or BP potential -----------------------------
1512 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1513 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
1515 C For the GB potential convert sigma'**2 into chi'
1518 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1522 write (iout,'(/a/)') 'Parameters of the BP potential:'
1523 write (iout,'(a/)') 'The epsilon array:'
1524 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1525 write (iout,'(/a)') 'One-body parameters:'
1526 write (iout,'(a,4x,5a)') 'residue',' sigma ','s||/s_|_^2',
1527 & ' chip0 ',' chip ',' alph '
1528 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),sigii(i),
1529 & chip0(i),chip(i),alp(i),i=1,ntyp)
1532 C--------------------- GBV potential -----------------------------------
1533 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1534 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1535 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1537 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1538 write (iout,'(a/)') 'The epsilon array:'
1539 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1540 write (iout,'(/a)') 'One-body parameters:'
1541 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1542 & 's||/s_|_^2',' chip ',' alph '
1543 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1544 & sigii(i),chip(i),alp(i),i=1,ntyp)
1547 C--------------------- Momo potential -----------------------------------
1549 read (isidep,*) (icharge(i),i=1,ntyp)
1550 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
1553 c! write (*,*) "Im in ", i, " ", j
1555 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
1556 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
1557 & chis(i,j),chis(j,i),
1558 & nstate(i,j),(wstate(k,i,j),k=1,4),
1563 & dtail(1,i,j),dtail(2,i,j),
1564 & epshead(i,j),sig0head(i,j),
1565 & rborn(i,j),rborn(j,i),
1566 & (wqdip(k,i,j),k=1,2),wquad(i,j),
1567 & alphapol(i,j),alphapol(j,i),
1568 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
1571 c! write (*,*) "nstate gly-gly", nstate(10,10)
1572 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
1573 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
1574 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
1577 c! DISABLE IT AT >>YOUR OWN PERIL<<
1582 sigma(i,j) = sigma(j,i)
1583 nstate(i,j) = nstate(j,i)
1584 sigmap1(i,j) = sigmap1(j,i)
1585 sigmap2(i,j) = sigmap2(j,i)
1586 sigiso1(i,j) = sigiso1(j,i)
1587 sigiso2(i,j) = sigiso2(j,i)
1590 alphasur(k,i,j) = alphasur(k,j,i)
1591 wstate(k,i,j) = wstate(k,j,i)
1592 alphiso(k,i,j) = alphiso(k,j,i)
1595 dhead(2,1,i,j) = dhead(1,1,j,i)
1596 dhead(2,2,i,j) = dhead(1,2,j,i)
1597 dhead(1,1,i,j) = dhead(2,1,j,i)
1598 dhead(1,2,i,j) = dhead(2,2,j,i)
1599 dtail(1,i,j) = dtail(1,j,i)
1600 dtail(2,i,j) = dtail(2,j,i)
1603 c! write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
1604 c! write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
1605 c! dhead(l,k,i,j) = dhead(k,l,j,i)
1609 epshead(i,j) = epshead(j,i)
1610 sig0head(i,j) = sig0head(j,i)
1613 wqdip(k,i,j) = wqdip(k,j,i)
1616 wquad(i,j) = wquad(j,i)
1617 epsintab(i,j) = epsintab(j,i)
1622 if (.not.lprint) goto 70
1624 & "Parameters of the new physics-based SC-SC interaction potential"
1625 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
1626 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
1629 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
1630 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
1631 & chipp(i,j),chipp(j,i)
1634 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
1635 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
1639 write (iout,'(2(a3,1x),8(0pf10.3))')
1640 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
1641 & sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
1644 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
1645 & ' wst2',' wst3',' wst4',' dh11',' dh21',
1646 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
1650 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
1651 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
1652 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1653 & epshead(i,j),sig0head(i,j)
1656 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
1657 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
1661 write (iout,'(2(a3,1x),2i4,5f10.3)')
1662 & restyp(i),restyp(j),icharge(i),icharge(j),
1663 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
1666 write (iout,'(/12a)') 'Residues',
1668 & ' alphpol2',' alphiso1',' alpiso2',
1669 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
1673 write (iout,'(2(a3,1x),11f10.3)')
1674 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
1675 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
1680 C For the GB potential convert sigma'**2 into chi'
1682 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1685 write (iout,'(/a/)') 'Parameters of the BP potential:'
1686 write (iout,'(a/)') 'The epsilon array:'
1687 CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1688 write (iout,'(/a)') 'One-body parameters:'
1689 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1691 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1692 & chip(i),alp(i),i=1,ntyp)
1696 C-----------------------------------------------------------------------
1697 C Calculate the "working" parameters of SC interactions.
1702 alphasur(k,j,i)=alphasur(k,i,j)
1703 alphiso(k,j,i)=alphiso(k,i,j)
1704 wstate(k,j,i)=wstate(k,i,j)
1707 wqdip(k,j,i)=wqdip(k,i,j)
1711 dhead(l,k,j,i)=dhead(l,k,i,j)
1719 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1720 sigma(j,i)=sigma(i,j)
1721 rs0(i,j)=dwa16*sigma(i,j)
1728 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1729 & 'Working parameters of the SC interactions:',
1730 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1735 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 )
1745 sigeps=dsign(1.0D0,epsij)
1747 aa(i,j)=epsij*rrij*rrij
1748 bb(i,j)=-sigeps*epsij*rrij
1751 IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
1752 sigt1sq=sigma0(i)**2
1753 sigt2sq=sigma0(j)**2
1756 ratsig1=sigt2sq/sigt1sq
1757 ratsig2=1.0D0/ratsig1
1758 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1759 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1760 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1764 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1765 sigmaii(i,j)=rsum_max
1766 sigmaii(j,i)=rsum_max
1768 c sigmaii(i,j)=r0(i,j)
1769 c sigmaii(j,i)=r0(i,j)
1771 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1772 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1773 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1774 augm(i,j)=epsij*r_augm**(2*expon)
1775 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1783 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1784 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1785 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1787 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1789 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1790 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1791 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1792 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1793 & chis(i,j),chis(j,i),
1794 & nstate(i,j),(wstate(k,i,j),k=1,4),
1795 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1796 & epshead(i,j),sig0head(i,j),
1797 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1798 & alphapol(i,j),alphapol(j,i),
1799 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1806 C Define the SC-p interaction constants
1810 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1812 c aad(i,1)=0.3D0*4.0D0**12
1813 C Following line for constants currently implemented
1814 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1815 aad(i,1)=1.5D0*4.0D0**12
1816 c aad(i,1)=0.17D0*5.6D0**12
1818 C "Soft" SC-p repulsion
1820 C Following line for constants currently implemented
1821 c aad(i,1)=0.3D0*4.0D0**6
1822 C "Hard" SC-p repulsion
1823 bad(i,1)=3.0D0*4.0D0**6
1824 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1833 C 8/9/01 Read the SC-p interaction constants from file
1836 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1839 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1840 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1841 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1842 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1846 write (iout,*) "Parameters of SC-p interactions:"
1848 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1849 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1854 C Define the constants of the disulfide bridge
1858 c Old arbitrary potential - commented out.
1863 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1864 c energy surface of diethyl disulfide.
1865 c A. Liwo and U. Kozlowska, 11/24/03
1874 write (iout,*) dyn_ss,'dyndyn'
1876 ss_depth=ebr/wsc-0.25*eps(1,1)
1877 C write(iout,*) akcm,whpb,wsc,'KURWA'
1878 Ht=Ht/wsc-0.25*eps(1,1)
1887 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1891 write (iout,'(/a)') "Disulfide bridge parameters:"
1892 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1893 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1894 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1895 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1898 if (shield_mode.gt.0) then
1900 C VSolvSphere the volume of solving sphere
1902 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1903 C there will be no distinction between proline peptide group and normal peptide
1904 C group in case of shielding parameters
1905 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1906 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1907 write (iout,*) VSolvSphere,VSolvSphere_div
1908 C long axis of side chain
1910 long_r_sidechain(i)=vbldsc0(1,i)
1911 short_r_sidechain(i)=sigma0(i)
1916 111 write (iout,*) "Error reading bending energy parameters."
1918 112 write (iout,*) "Error reading rotamer energy parameters."
1920 113 write (iout,*) "Error reading torsional energy parameters."
1922 114 write (iout,*) "Error reading double torsional energy parameters."
1925 & "Error reading cumulant (multibody energy) parameters."
1927 116 write (iout,*) "Error reading electrostatic energy parameters."
1929 1161 write (iout,*) "Error reading lipid energy parameters."
1931 117 write (iout,*) "Error reading side chain interaction parameters."
1933 118 write (iout,*) "Error reading SCp interaction parameters."
1935 119 write (iout,*) "Error reading SCCOR parameters"
1937 121 write (iout,*) "Error reading bond parameters"
1940 call MPI_Finalize(Ierror)