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)
681 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
682 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
683 itype2loc(ntyp1)=nloctyp
684 iloctyp(nloctyp)=ntyp1
686 itype2loc(-i)=-itype2loc(i)
695 iloctyp(-i)=-iloctyp(i)
697 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
698 c write (iout,*) "nloctyp",nloctyp,
699 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
702 c write (iout,*) "NEWCORR",i
703 read (ifourier,*,end=115,err=115)
706 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
709 c write (iout,*) "NEWCORR BNEW1"
710 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
713 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
716 c write (iout,*) "NEWCORR BNEW2"
717 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
719 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
720 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
722 c write (iout,*) "NEWCORR CCNEW"
723 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
725 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
726 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
728 c write (iout,*) "NEWCORR DDNEW"
729 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
733 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
737 c write (iout,*) "NEWCORR EENEW1"
738 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
740 read (ifourier,*,end=115,err=115) e0new(ii,i)
742 c write (iout,*) (e0new(ii,i),ii=1,3)
744 c write (iout,*) "NEWCORR EENEW"
747 ccnew(ii,1,i)=ccnew(ii,1,i)/2
748 ccnew(ii,2,i)=ccnew(ii,2,i)/2
749 ddnew(ii,1,i)=ddnew(ii,1,i)/2
750 ddnew(ii,2,i)=ddnew(ii,2,i)/2
755 bnew1(ii,1,-i)= bnew1(ii,1,i)
756 bnew1(ii,2,-i)=-bnew1(ii,2,i)
757 bnew2(ii,1,-i)= bnew2(ii,1,i)
758 bnew2(ii,2,-i)=-bnew2(ii,2,i)
761 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
762 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
763 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
764 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
765 ccnew(ii,1,-i)=ccnew(ii,1,i)
766 ccnew(ii,2,-i)=-ccnew(ii,2,i)
767 ddnew(ii,1,-i)=ddnew(ii,1,i)
768 ddnew(ii,2,-i)=-ddnew(ii,2,i)
770 e0new(1,-i)= e0new(1,i)
771 e0new(2,-i)=-e0new(2,i)
772 e0new(3,-i)=-e0new(3,i)
774 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
775 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
776 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
777 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
781 write (iout,'(a)') "Coefficients of the multibody terms"
782 do i=-nloctyp+1,nloctyp-1
783 write (iout,*) "Type: ",onelet(iloctyp(i))
784 write (iout,*) "Coefficients of the expansion of B1"
786 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
788 write (iout,*) "Coefficients of the expansion of B2"
790 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
792 write (iout,*) "Coefficients of the expansion of C"
793 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
794 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
795 write (iout,*) "Coefficients of the expansion of D"
796 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
797 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
798 write (iout,*) "Coefficients of the expansion of E"
799 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
802 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
807 IF (SPLIT_FOURIERTOR) THEN
809 c write (iout,*) "NEWCORR TOR",i
810 read (ifourier,*,end=115,err=115)
813 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
816 c write (iout,*) "NEWCORR BNEW1 TOR"
817 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
820 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
823 c write (iout,*) "NEWCORR BNEW2 TOR"
824 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
826 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
827 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
829 c write (iout,*) "NEWCORR CCNEW TOR"
830 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
832 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
833 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
835 c write (iout,*) "NEWCORR DDNEW TOR"
836 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
840 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
844 c write (iout,*) "NEWCORR EENEW1 TOR"
845 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
847 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
849 c write (iout,*) (e0newtor(ii,i),ii=1,3)
851 c write (iout,*) "NEWCORR EENEW TOR"
854 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
855 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
856 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
857 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
862 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
863 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
864 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
865 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
868 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
869 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
870 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
871 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
873 e0newtor(1,-i)= e0newtor(1,i)
874 e0newtor(2,-i)=-e0newtor(2,i)
875 e0newtor(3,-i)=-e0newtor(3,i)
877 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
878 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
879 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
880 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
885 & "Single-body coefficients of the torsional potentials"
886 do i=-nloctyp+1,nloctyp-1
887 write (iout,*) "Type: ",onelet(iloctyp(i))
888 write (iout,*) "Coefficients of the expansion of B1tor"
890 write (iout,'(3hB1(,i1,1h),3f10.5)')
891 & j,(bnew1tor(k,j,i),k=1,3)
893 write (iout,*) "Coefficients of the expansion of B2tor"
895 write (iout,'(3hB2(,i1,1h),3f10.5)')
896 & j,(bnew2tor(k,j,i),k=1,3)
898 write (iout,*) "Coefficients of the expansion of Ctor"
899 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
900 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
901 write (iout,*) "Coefficients of the expansion of Dtor"
902 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
903 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
904 write (iout,*) "Coefficients of the expansion of Etor"
905 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
908 write (iout,'(1hE,2i1,2f10.5)')
909 & j,k,(eenewtor(l,j,k,i),l=1,2)
915 do i=-nloctyp+1,nloctyp-1
918 bnew1tor(ii,j,i)=bnew1(ii,j,i)
923 bnew2tor(ii,j,i)=bnew2(ii,j,i)
927 ccnewtor(ii,1,i)=ccnew(ii,1,i)
928 ccnewtor(ii,2,i)=ccnew(ii,2,i)
929 ddnewtor(ii,1,i)=ddnew(ii,1,i)
930 ddnewtor(ii,2,i)=ddnew(ii,2,i)
933 ENDIF !SPLIT_FOURIER_TOR
936 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
938 read (ifourier,*,end=115,err=115)
939 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
941 write (iout,*) 'Type ',onelet(iloctyp(i))
942 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
956 c B1tilde(1,i) = b(3)
957 c B1tilde(2,i) =-b(5)
958 c B1tilde(1,-i) =-b(3)
959 c B1tilde(2,-i) =b(5)
966 cc B1tilde(1,i) = b(3,i)
967 cc B1tilde(2,i) =-b(5,i)
968 C B1tilde(1,-i) =-b(3,i)
969 C B1tilde(2,-i) =b(5,i)
970 cc b1tilde(1,i)=0.0d0
971 cc b1tilde(2,i)=0.0d0
983 CCold(1,1,-i)= b(7,i)
984 CCold(2,2,-i)=-b(7,i)
985 CCold(2,1,-i)=-b(9,i)
986 CCold(1,2,-i)=-b(9,i)
991 c Ctilde(1,1,i)= CCold(1,1,i)
992 c Ctilde(1,2,i)= CCold(1,2,i)
993 c Ctilde(2,1,i)=-CCold(2,1,i)
994 c Ctilde(2,2,i)=-CCold(2,2,i)
996 c Ctilde(1,1,i)=0.0d0
997 c Ctilde(1,2,i)=0.0d0
998 c Ctilde(2,1,i)=0.0d0
999 c Ctilde(2,2,i)=0.0d0
1000 DDold(1,1,i)= b(6,i)
1001 DDold(2,2,i)=-b(6,i)
1002 DDold(2,1,i)= b(8,i)
1003 DDold(1,2,i)= b(8,i)
1004 DDold(1,1,-i)= b(6,i)
1005 DDold(2,2,-i)=-b(6,i)
1006 DDold(2,1,-i)=-b(8,i)
1007 DDold(1,2,-i)=-b(8,i)
1012 c Dtilde(1,1,i)= DD(1,1,i)
1013 c Dtilde(1,2,i)= DD(1,2,i)
1014 c Dtilde(2,1,i)=-DD(2,1,i)
1015 c Dtilde(2,2,i)=-DD(2,2,i)
1017 c Dtilde(1,1,i)=0.0d0
1018 c Dtilde(1,2,i)=0.0d0
1019 c Dtilde(2,1,i)=0.0d0
1020 c Dtilde(2,2,i)=0.0d0
1021 EEold(1,1,i)= b(10,i)+b(11,i)
1022 EEold(2,2,i)=-b(10,i)+b(11,i)
1023 EEold(2,1,i)= b(12,i)-b(13,i)
1024 EEold(1,2,i)= b(12,i)+b(13,i)
1025 EEold(1,1,-i)= b(10,i)+b(11,i)
1026 EEold(2,2,-i)=-b(10,i)+b(11,i)
1027 EEold(2,1,-i)=-b(12,i)+b(13,i)
1028 EEold(1,2,-i)=-b(12,i)-b(13,i)
1029 write(iout,*) "TU DOCHODZE"
1035 c ee(2,1,i)=ee(1,2,i)
1040 &"Coefficients of the cumulants (independent of valence angles)"
1041 do i=-nloctyp+1,nloctyp-1
1042 write (iout,*) 'Type ',onelet(iloctyp(i))
1044 write(iout,'(2f10.5)') B(3,i),B(5,i)
1046 write(iout,'(2f10.5)') B(2,i),B(4,i)
1049 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1053 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1057 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1062 C write (iout,*) 'KURWAKURWA'
1065 C Read torsional parameters in old format
1067 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1068 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1069 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1072 read (itorp,'(a)',end=113,err=113)
1074 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1080 write (iout,'(/a/)') 'Torsional constants:'
1083 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1084 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1090 C Read torsional parameters
1092 IF (TOR_MODE.eq.0) THEN
1094 read (itorp,*,end=113,err=113) ntortyp
1095 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1096 write (iout,*) 'ntortyp',ntortyp
1098 itype2loc(i)=itortyp(i)
1101 itype2loc(-i)=-itype2loc(i)
1105 itortyp(i)=-itortyp(-i)
1107 c write (iout,*) 'ntortyp',ntortyp
1109 do j=-ntortyp+1,ntortyp-1
1110 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1112 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1113 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1116 do k=1,nterm(i,j,iblock)
1117 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1119 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1120 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1121 v0ij=v0ij+si*v1(k,i,j,iblock)
1124 do k=1,nlor(i,j,iblock)
1125 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1126 & vlor2(k,i,j),vlor3(k,i,j)
1127 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1130 v0(-i,-j,iblock)=v0ij
1136 write (iout,'(/a/)') 'Torsional constants:'
1140 write (iout,*) 'ityp',i,' jtyp',j
1141 write (iout,*) 'Fourier constants'
1142 do k=1,nterm(i,j,iblock)
1143 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1146 write (iout,*) 'Lorenz constants'
1147 do k=1,nlor(i,j,iblock)
1148 write (iout,'(3(1pe15.5))')
1149 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1156 C 6/23/01 Read parameters for double torsionals
1160 do j=-ntortyp+1,ntortyp-1
1161 do k=-ntortyp+1,ntortyp-1
1162 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1163 c write (iout,*) "OK onelett",
1166 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1167 & .or. t3.ne.toronelet(k)) then
1168 write (iout,*) "Error in double torsional parameter file",
1171 call MPI_Finalize(Ierror)
1173 stop "Error in double torsional parameter file"
1175 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1176 & ntermd_2(i,j,k,iblock)
1177 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1178 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1179 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1180 & ntermd_1(i,j,k,iblock))
1181 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1182 & ntermd_1(i,j,k,iblock))
1183 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1184 & ntermd_1(i,j,k,iblock))
1185 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1186 & ntermd_1(i,j,k,iblock))
1187 C Martix of D parameters for one dimesional foureir series
1188 do l=1,ntermd_1(i,j,k,iblock)
1189 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1190 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1191 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1192 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1193 c write(iout,*) "whcodze" ,
1194 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1196 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1197 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1198 & v2s(m,l,i,j,k,iblock),
1199 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1200 C Martix of D parameters for two dimesional fourier series
1201 do l=1,ntermd_2(i,j,k,iblock)
1203 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1204 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1205 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1206 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1215 write (iout,*) 'Constants for double torsionals'
1218 do j=-ntortyp+1,ntortyp-1
1219 do k=-ntortyp+1,ntortyp-1
1220 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1221 & ' nsingle',ntermd_1(i,j,k,iblock),
1222 & ' ndouble',ntermd_2(i,j,k,iblock)
1224 write (iout,*) 'Single angles:'
1225 do l=1,ntermd_1(i,j,k,iblock)
1226 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1227 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1228 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1229 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1232 write (iout,*) 'Pairs of angles:'
1233 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1234 do l=1,ntermd_2(i,j,k,iblock)
1235 write (iout,'(i5,20f10.5)')
1236 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1239 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1240 do l=1,ntermd_2(i,j,k,iblock)
1241 write (iout,'(i5,20f10.5)')
1242 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1243 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1252 ELSE IF (TOR_MODE.eq.1) THEN
1254 C read valence-torsional parameters
1255 read (itorp,*,end=113,err=113) ntortyp
1257 write (iout,*) "Valence-torsional parameters read in ntortyp",
1259 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1260 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1262 itortyp(i)=-itortyp(-i)
1264 do i=-ntortyp+1,ntortyp-1
1265 do j=-ntortyp+1,ntortyp-1
1266 C first we read the cos and sin gamma parameters
1267 read (itorp,'(13x,a)',end=113,err=113) string
1268 write (iout,*) i,j,string
1269 read (itorp,*,end=113,err=113)
1270 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1271 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1272 do k=1,nterm_kcc(j,i)
1273 do l=1,nterm_kcc_Tb(j,i)
1274 do ll=1,nterm_kcc_Tb(j,i)
1275 read (itorp,*,end=113,err=113) ii,jj,kk,
1276 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1284 c AL 4/8/16: Calculate coefficients from one-body parameters
1287 itortyp(i)=itype2loc(i)
1291 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1293 do i=-ntortyp+1,ntortyp-1
1294 do j=-ntortyp+1,ntortyp-1
1297 do k=1,nterm_kcc_Tb(j,i)
1298 do l=1,nterm_kcc_Tb(j,i)
1299 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1300 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1301 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1302 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1305 do k=1,nterm_kcc_Tb(j,i)
1306 do l=1,nterm_kcc_Tb(j,i)
1307 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1308 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1309 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1310 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1313 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)
1319 if (tor_mode.gt.0 .and. lprint) then
1320 c Print valence-torsional parameters
1322 & "Parameters of the valence-torsional potentials"
1323 do i=-ntortyp+1,ntortyp-1
1324 do j=-ntortyp+1,ntortyp-1
1325 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1326 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1327 do k=1,nterm_kcc(j,i)
1328 do l=1,nterm_kcc_Tb(j,i)
1329 do ll=1,nterm_kcc_Tb(j,i)
1330 write (iout,'(3i5,2f15.4)')
1331 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1340 C Read of Side-chain backbone correlation parameters
1341 C Modified 11 May 2012 by Adasko
1344 read (isccor,*,end=119,err=119) nsccortyp
1345 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1347 isccortyp(i)=-isccortyp(-i)
1349 iscprol=isccortyp(20)
1350 c write (iout,*) 'ntortyp',ntortyp
1352 cc maxinter is maximum interaction sites
1356 read (isccor,*,end=119,err=119)
1357 &nterm_sccor(i,j),nlor_sccor(i,j)
1358 c write (iout,*) nterm_sccor(i,j)
1364 nterm_sccor(-i,j)=nterm_sccor(i,j)
1365 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1366 nterm_sccor(i,-j)=nterm_sccor(i,j)
1367 c write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1368 c & nterm_sccor(-i,-j),nterm_sccor(i,-j)
1369 do k=1,nterm_sccor(i,j)
1370 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1372 if (j.eq.iscprol) then
1373 if (i.eq.isccortyp(10)) then
1374 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1375 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1377 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1378 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1379 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1380 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1381 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1382 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1383 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1384 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1387 if (i.eq.isccortyp(10)) then
1388 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1389 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1391 if (j.eq.isccortyp(10)) then
1392 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1393 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1395 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1396 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1397 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1398 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1399 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1400 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1404 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1405 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1406 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1407 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1410 do k=1,nlor_sccor(i,j)
1411 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1412 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1413 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1414 &(1+vlor3sccor(k,i,j)**2)
1416 v0sccor(l,i,j)=v0ijsccor
1417 v0sccor(l,-i,j)=v0ijsccor1
1418 v0sccor(l,i,-j)=v0ijsccor2
1419 v0sccor(l,-i,-j)=v0ijsccor3
1425 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1429 write (iout,*) 'ityp',i,' jtyp',j
1430 write (iout,*) 'Fourier constants'
1431 do k=1,nterm_sccor(i,j)
1432 write (iout,'(2(1pe15.5))')
1433 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1435 write (iout,*) 'Lorenz constants'
1436 do k=1,nlor_sccor(i,j)
1437 write (iout,'(3(1pe15.5))')
1438 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1445 C Read electrostatic-interaction parameters
1448 write (iout,'(/a)') 'Electrostatic interaction constants:'
1449 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1450 & 'IT','JT','APP','BPP','AEL6','AEL3'
1452 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1453 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1454 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1455 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1460 app (i,j)=epp(i,j)*rri*rri
1461 bpp (i,j)=-2.0D0*epp(i,j)*rri
1462 ael6(i,j)=elpp6(i,j)*4.2D0**6
1463 ael3(i,j)=elpp3(i,j)*4.2D0**3
1465 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1466 & ael6(i,j),ael3(i,j)
1470 C Read side-chain interaction parameters.
1472 read (isidep,*) ipot,expon
1473 if (ipot.lt.1 .or. ipot.gt.6) then
1474 write (iout,'(2a)') 'Error while reading SC interaction',
1475 & 'potential file - unknown potential type.'
1479 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1480 & ', exponents are ',expon,2*expon
1481 goto (10,20,30,30,40,50,60) ipot
1482 C----------------------- LJ potential ---------------------------------
1483 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1485 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1486 write (iout,'(a/)') 'The epsilon array:'
1487 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1488 write (iout,'(/a)') 'One-body parameters:'
1489 write (iout,'(a,4x,a)') 'residue','sigma'
1490 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1493 C----------------------- LJK potential --------------------------------
1494 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1495 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1497 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1498 write (iout,'(a/)') 'The epsilon array:'
1499 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1500 write (iout,'(/a)') 'One-body parameters:'
1501 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1502 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1506 C---------------------- GB or BP potential -----------------------------
1507 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1508 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
1510 C For the GB potential convert sigma'**2 into chi'
1513 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1517 write (iout,'(/a/)') 'Parameters of the BP potential:'
1518 write (iout,'(a/)') 'The epsilon array:'
1519 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1520 write (iout,'(/a)') 'One-body parameters:'
1521 write (iout,'(a,4x,5a)') 'residue',' sigma ','s||/s_|_^2',
1522 & ' chip0 ',' chip ',' alph '
1523 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),sigii(i),
1524 & chip0(i),chip(i),alp(i),i=1,ntyp)
1527 C--------------------- GBV potential -----------------------------------
1528 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1529 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1530 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1532 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1533 write (iout,'(a/)') 'The epsilon array:'
1534 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1535 write (iout,'(/a)') 'One-body parameters:'
1536 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1537 & 's||/s_|_^2',' chip ',' alph '
1538 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1539 & sigii(i),chip(i),alp(i),i=1,ntyp)
1542 C--------------------- Momo potential -----------------------------------
1544 read (isidep,*) (icharge(i),i=1,ntyp)
1545 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
1548 c! write (*,*) "Im in ", i, " ", j
1550 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
1551 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
1552 & chis(i,j),chis(j,i),
1553 & nstate(i,j),(wstate(k,i,j),k=1,4),
1558 & dtail(1,i,j),dtail(2,i,j),
1559 & epshead(i,j),sig0head(i,j),
1560 & rborn(i,j),rborn(j,i),
1561 & (wqdip(k,i,j),k=1,2),wquad(i,j),
1562 & alphapol(i,j),alphapol(j,i),
1563 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
1566 c! write (*,*) "nstate gly-gly", nstate(10,10)
1567 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
1568 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
1569 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
1572 c! DISABLE IT AT >>YOUR OWN PERIL<<
1577 sigma(i,j) = sigma(j,i)
1578 nstate(i,j) = nstate(j,i)
1579 sigmap1(i,j) = sigmap1(j,i)
1580 sigmap2(i,j) = sigmap2(j,i)
1581 sigiso1(i,j) = sigiso1(j,i)
1582 sigiso2(i,j) = sigiso2(j,i)
1585 alphasur(k,i,j) = alphasur(k,j,i)
1586 wstate(k,i,j) = wstate(k,j,i)
1587 alphiso(k,i,j) = alphiso(k,j,i)
1590 dhead(2,1,i,j) = dhead(1,1,j,i)
1591 dhead(2,2,i,j) = dhead(1,2,j,i)
1592 dhead(1,1,i,j) = dhead(2,1,j,i)
1593 dhead(1,2,i,j) = dhead(2,2,j,i)
1594 dtail(1,i,j) = dtail(1,j,i)
1595 dtail(2,i,j) = dtail(2,j,i)
1598 c! write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
1599 c! write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
1600 c! dhead(l,k,i,j) = dhead(k,l,j,i)
1604 epshead(i,j) = epshead(j,i)
1605 sig0head(i,j) = sig0head(j,i)
1608 wqdip(k,i,j) = wqdip(k,j,i)
1611 wquad(i,j) = wquad(j,i)
1612 epsintab(i,j) = epsintab(j,i)
1617 if (.not.lprint) goto 70
1619 & "Parameters of the new physics-based SC-SC interaction potential"
1620 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
1621 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
1624 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
1625 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
1626 & chipp(i,j),chipp(j,i)
1629 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
1630 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
1634 write (iout,'(2(a3,1x),8(0pf10.3))')
1635 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
1636 & sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
1639 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
1640 & ' wst2',' wst3',' wst4',' dh11',' dh21',
1641 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
1645 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
1646 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
1647 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1648 & epshead(i,j),sig0head(i,j)
1651 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
1652 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
1656 write (iout,'(2(a3,1x),2i4,5f10.3)')
1657 & restyp(i),restyp(j),icharge(i),icharge(j),
1658 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
1661 write (iout,'(/12a)') 'Residues',
1663 & ' alphpol2',' alphiso1',' alpiso2',
1664 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
1668 write (iout,'(2(a3,1x),11f10.3)')
1669 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
1670 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
1675 C For the GB potential convert sigma'**2 into chi'
1677 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1680 write (iout,'(/a/)') 'Parameters of the BP potential:'
1681 write (iout,'(a/)') 'The epsilon array:'
1682 CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1683 write (iout,'(/a)') 'One-body parameters:'
1684 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1686 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1687 & chip(i),alp(i),i=1,ntyp)
1691 C-----------------------------------------------------------------------
1692 C Calculate the "working" parameters of SC interactions.
1697 alphasur(k,j,i)=alphasur(k,i,j)
1698 alphiso(k,j,i)=alphiso(k,i,j)
1699 wstate(k,j,i)=wstate(k,i,j)
1702 wqdip(k,j,i)=wqdip(k,i,j)
1706 dhead(l,k,j,i)=dhead(l,k,i,j)
1714 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1715 sigma(j,i)=sigma(i,j)
1716 rs0(i,j)=dwa16*sigma(i,j)
1723 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1724 & 'Working parameters of the SC interactions:',
1725 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1730 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 )
1740 sigeps=dsign(1.0D0,epsij)
1742 aa(i,j)=epsij*rrij*rrij
1743 bb(i,j)=-sigeps*epsij*rrij
1746 IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
1747 sigt1sq=sigma0(i)**2
1748 sigt2sq=sigma0(j)**2
1751 ratsig1=sigt2sq/sigt1sq
1752 ratsig2=1.0D0/ratsig1
1753 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1754 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1755 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1759 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1760 sigmaii(i,j)=rsum_max
1761 sigmaii(j,i)=rsum_max
1763 c sigmaii(i,j)=r0(i,j)
1764 c sigmaii(j,i)=r0(i,j)
1766 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1767 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1768 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1769 augm(i,j)=epsij*r_augm**(2*expon)
1770 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1778 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1779 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1780 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1782 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
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),
1786 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1787 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1788 & chis(i,j),chis(j,i),
1789 & nstate(i,j),(wstate(k,i,j),k=1,4),
1790 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1791 & epshead(i,j),sig0head(i,j),
1792 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1793 & alphapol(i,j),alphapol(j,i),
1794 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1801 C Define the SC-p interaction constants
1805 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1807 c aad(i,1)=0.3D0*4.0D0**12
1808 C Following line for constants currently implemented
1809 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1810 aad(i,1)=1.5D0*4.0D0**12
1811 c aad(i,1)=0.17D0*5.6D0**12
1813 C "Soft" SC-p repulsion
1815 C Following line for constants currently implemented
1816 c aad(i,1)=0.3D0*4.0D0**6
1817 C "Hard" SC-p repulsion
1818 bad(i,1)=3.0D0*4.0D0**6
1819 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1828 C 8/9/01 Read the SC-p interaction constants from file
1831 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1834 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1835 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1836 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1837 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1841 write (iout,*) "Parameters of SC-p interactions:"
1843 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1844 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1849 C Define the constants of the disulfide bridge
1853 c Old arbitrary potential - commented out.
1858 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1859 c energy surface of diethyl disulfide.
1860 c A. Liwo and U. Kozlowska, 11/24/03
1869 write (iout,*) dyn_ss,'dyndyn'
1871 ss_depth=ebr/wsc-0.25*eps(1,1)
1872 C write(iout,*) akcm,whpb,wsc,'KURWA'
1873 Ht=Ht/wsc-0.25*eps(1,1)
1882 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1886 write (iout,'(/a)') "Disulfide bridge parameters:"
1887 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1888 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1889 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1890 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1893 if (shield_mode.gt.0) then
1895 C VSolvSphere the volume of solving sphere
1897 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1898 C there will be no distinction between proline peptide group and normal peptide
1899 C group in case of shielding parameters
1900 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1901 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1902 write (iout,*) VSolvSphere,VSolvSphere_div
1903 C long axis of side chain
1905 long_r_sidechain(i)=vbldsc0(1,i)
1906 short_r_sidechain(i)=sigma0(i)
1911 111 write (iout,*) "Error reading bending energy parameters."
1913 112 write (iout,*) "Error reading rotamer energy parameters."
1915 113 write (iout,*) "Error reading torsional energy parameters."
1917 114 write (iout,*) "Error reading double torsional energy parameters."
1920 & "Error reading cumulant (multibody energy) parameters."
1922 116 write (iout,*) "Error reading electrostatic energy parameters."
1924 1161 write (iout,*) "Error reading lipid energy parameters."
1926 117 write (iout,*) "Error reading side chain interaction parameters."
1928 118 write (iout,*) "Error reading SCp interaction parameters."
1930 119 write (iout,*) "Error reading SCCOR parameters"
1932 121 write (iout,*) "Error reading bond parameters"
1935 call MPI_Finalize(Ierror)