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'
26 include 'COMMON.LANGEVIN'
28 character*1 onelett(4) /"G","A","P","D"/
29 character*1 toronelet(-2:2) /"p","a","G","A","P"/
31 dimension blower(3,3,maxlob)
32 character*800 controlcard
33 character*256 bondname_t,thetname_t,rotname_t,torname_t,
34 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
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),ww0(i))
60 write (iout,*) i,key(:ilen(key)),ww(i)
62 call reada(controlcard,'SCAL14',scal14,0.4d0)
63 call reada(controlcard,'SCALSCP',scalscp,1.0d0)
64 call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
65 call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
66 call reada(controlcard,'LIPSCALE',lipscale,1.0D0)
67 c write (iout,*) "lipscale",lipscale
68 r0_corr=cutoff_corr-delt_corr
70 write (iout,*) "iparm",iparm," myparm",myparm
71 c If reading not own parameters, skip assignment
72 call reada(controlcard,"D0CM",d0cm,3.78d0)
73 call reada(controlcard,"AKCM",akcm,15.1d0)
74 call reada(controlcard,"AKTH",akth,11.0d0)
75 call reada(controlcard,"AKCT",akct,12.0d0)
76 call reada(controlcard,"V1SS",v1ss,-1.08d0)
77 call reada(controlcard,"V2SS",v2ss,7.61d0)
78 call reada(controlcard,"V3SS",v3ss,13.7d0)
79 call reada(controlcard,"EBR",ebr,-5.50D0)
80 call reada(controlcard,"DTRISS",dtriss,1.0D0)
81 call reada(controlcard,"ATRISS",atriss,0.3D0)
82 call reada(controlcard,"BTRISS",btriss,0.02D0)
83 call reada(controlcard,"CTRISS",ctriss,1.0D0)
84 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
86 C dyn_ss_mask(i)=.false.
90 c Old arbitrary potential - commented out.
95 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
96 c energy surface of diethyl disulfide.
97 c A. Liwo and U. Kozlowska, 11/24/03
109 dyn_ssbond_ij(i,j)=1.0d300
112 call reada(controlcard,"HT",Ht,0.0D0)
114 C ss_depth=ebr/wsc-0.25*eps(1,1)
115 C write(iout,*) HT,wsc,eps(1,1),'KURWA'
116 C Ht=Ht/wsc-0.25*eps(1,1)
125 C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
128 if (iparm.eq.myparm .or. .not.separate_parset) then
131 c Setup weights for UNRES
155 c write (iout,*) "PARMREAD: wsaxs",wsaxs
160 c write(iout,*)"PARMREAD: wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
161 c & " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
165 C Juyong:READ init_vars
166 C Initialize variables!
167 C Juyong:READ read_info
168 C READ fragment information!!
169 C both routines should be in dfa.F file!!
170 write (iout,*) "Before initializing DFA",wdfa_dist,wdfa_tor,
172 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
173 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
174 write (iout,*) "Calling init_dfa_vars"
177 write (iout,*) 'init_dfa_vars finished!'
180 write (iout,*) 'read_dfa_info finished!'
184 call card_concat(controlcard,.false.)
186 c Return if not own parameters
188 if (iparm.ne.myparm .and. separate_parset) return
190 call reads(controlcard,"BONDPAR",bondname_t,bondname)
191 open (ibond,file=bondname_t,status='old')
193 call reads(controlcard,"THETPAR",thetname_t,thetname)
194 open (ithep,file=thetname_t,status='old')
196 call reads(controlcard,"ROTPAR",rotname_t,rotname)
197 open (irotam,file=rotname_t,status='old')
199 call reads(controlcard,"TORPAR",torname_t,torname)
200 open (itorp,file=torname_t,status='old')
202 call reads(controlcard,"TORDPAR",tordname_t,tordname)
203 write (iout,*) "tor_mode",tor_mode
206 & open (itordp,file=tordname_t,status='old')
208 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
209 open (isccor,file=sccorname_t,status='old')
211 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
212 open (ifourier,file=fouriername_t,status='old')
214 call reads(controlcard,"ELEPAR",elename_t,elename)
215 open (ielep,file=elename_t,status='old')
217 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
218 open (isidep,file=sidename_t,status='old')
220 call reads(controlcard,"SCPPAR",scpname_t,scpname)
221 open (iscpp,file=scpname_t,status='old')
223 write (iout,*) "Parameter set:",iparm
224 write (iout,*) "Energy-term weights:"
226 write (iout,'(a16,f10.5)') wname(i),ww(i)
228 write (iout,*) "Sidechain potential file : ",
229 & sidename_t(:ilen(sidename_t))
231 write (iout,*) "SCp potential file : ",
232 & scpname_t(:ilen(scpname_t))
234 write (iout,*) "Electrostatic potential file : ",
235 & elename_t(:ilen(elename_t))
236 write (iout,*) "Cumulant coefficient file : ",
237 & fouriername_t(:ilen(fouriername_t))
238 write (iout,*) "Torsional parameter file : ",
239 & torname_t(:ilen(torname_t))
240 write (iout,*) "Double torsional parameter file : ",
241 & tordname_t(:ilen(tordname_t))
242 write (iout,*) "Backbone-rotamer parameter file : ",
243 & sccorname(:ilen(sccorname))
244 write (iout,*) "Bond & inertia constant file : ",
245 & bondname_t(:ilen(bondname_t))
246 write (iout,*) "Bending parameter file : ",
247 & thetname_t(:ilen(thetname_t))
248 write (iout,*) "Rotamer parameter file : ",
249 & rotname_t(:ilen(rotname_t))
252 c Read the virtual-bond parameters, masses, and moments of inertia
253 c and Stokes' radii of the peptide group and side chains
256 read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp,mp,ip,pstok
259 read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i),
260 & msc(i),isc(i),restok(i)
261 dsc(i) = vbldsc0(1,i)
265 dsc_inv(i)=1.0D0/dsc(i)
269 read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk,
272 read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
273 & aksc(j,i),abond0(j,i),j=1,nbondterm(i)),msc(i),isc(i),restok(i)
274 dsc(i) = vbldsc0(1,i)
278 dsc_inv(i)=1.0D0/dsc(i)
283 write(iout,'(/a/)')"Force constants virtual bonds:"
284 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
286 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
288 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
289 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
291 write (iout,'(13x,3f10.5)')
292 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
296 c write (iout,*) "iliptranpar",iliptranpar
297 c write (iout,*) "liptranname ",liptranname
298 read(iliptranpar,*,end=1161,err=1161) pepliptran
300 read(iliptranpar,*,end=1161,err=1161) liptranene(i)
304 write (iout,'(/a)') "Water-lipid transfer parameters"
305 write (iout,'(a3,3x,f10.5)') 'p',pepliptran
307 write (iout,'(a3,3x,f10.5)') restyp(i),liptranene(i)
312 C Read the parameters of the probability distribution/energy expression
313 C of the virtual-bond valence angles theta
316 read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
317 & (bthet(j,i,1,1),j=1,2)
318 read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
319 read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
320 read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
324 athet(1,i,1,-1)=athet(1,i,1,1)
325 athet(2,i,1,-1)=athet(2,i,1,1)
326 bthet(1,i,1,-1)=-bthet(1,i,1,1)
327 bthet(2,i,1,-1)=-bthet(2,i,1,1)
328 athet(1,i,-1,1)=-athet(1,i,1,1)
329 athet(2,i,-1,1)=-athet(2,i,1,1)
330 bthet(1,i,-1,1)=bthet(1,i,1,1)
331 bthet(2,i,-1,1)=bthet(2,i,1,1)
335 athet(1,i,-1,-1)=athet(1,-i,1,1)
336 athet(2,i,-1,-1)=-athet(2,-i,1,1)
337 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
338 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
339 athet(1,i,-1,1)=athet(1,-i,1,1)
340 athet(2,i,-1,1)=-athet(2,-i,1,1)
341 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
342 bthet(2,i,-1,1)=bthet(2,-i,1,1)
343 athet(1,i,1,-1)=-athet(1,-i,1,1)
344 athet(2,i,1,-1)=athet(2,-i,1,1)
345 bthet(1,i,1,-1)=bthet(1,-i,1,1)
346 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
351 polthet(j,i)=polthet(j,-i)
354 gthet(j,i)=gthet(j,-i)
360 c & 'Parameters of the virtual-bond valence angles:'
361 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
362 c & ' ATHETA0 ',' A1 ',' A2 ',
365 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
366 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
368 c write (iout,'(/a/9x,5a/79(1h-))')
369 c & 'Parameters of the expression for sigma(theta_c):',
370 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
371 c & ' ALPH3 ',' SIGMA0C '
373 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
374 c & (polthet(j,i),j=0,3),sigc0(i)
376 c write (iout,'(/a/9x,5a/79(1h-))')
377 c & 'Parameters of the second gaussian:',
378 c & ' THETA0 ',' SIGMA0 ',' G1 ',
381 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
382 c & sig0(i),(gthet(j,i),j=1,3)
385 & 'Parameters of the virtual-bond valence angles:'
386 write (iout,'(/a/9x,5a/79(1h-))')
387 & 'Coefficients of expansion',
388 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
389 & ' b1*10^1 ',' b2*10^1 '
391 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
392 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
393 & (10*bthet(j,i,1,1),j=1,2)
395 write (iout,'(/a/9x,5a/79(1h-))')
396 & 'Parameters of the expression for sigma(theta_c):',
397 & ' alpha0 ',' alph1 ',' alph2 ',
398 & ' alhp3 ',' sigma0c '
400 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
401 & (polthet(j,i),j=0,3),sigc0(i)
403 write (iout,'(/a/9x,5a/79(1h-))')
404 & 'Parameters of the second gaussian:',
405 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
408 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
409 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
413 IF (tor_mode.eq.0) THEN
415 C Read the parameters of Utheta determined from ab initio surfaces
416 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
418 read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
419 & ntheterm3,nsingle,ndouble
420 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
421 read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
423 ithetyp(i)=-ithetyp(-i)
426 do i=-maxthetyp,maxthetyp
427 do j=-maxthetyp,maxthetyp
428 do k=-maxthetyp,maxthetyp
429 aa0thet(i,j,k,iblock)=0.0d0
431 aathet(l,i,j,k,iblock)=0.0d0
435 bbthet(m,l,i,j,k,iblock)=0.0d0
436 ccthet(m,l,i,j,k,iblock)=0.0d0
437 ddthet(m,l,i,j,k,iblock)=0.0d0
438 eethet(m,l,i,j,k,iblock)=0.0d0
444 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
445 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
453 C write (iout,*) "KURWA1"
456 do j=-nthetyp,nthetyp
457 do k=-nthetyp,nthetyp
458 read (ithep,'(6a)',end=111,err=111) res1
459 write(iout,*) res1,i,j,k
460 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
461 read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
463 read (ithep,*,end=111,err=111)
464 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
465 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
466 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
467 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
469 read (ithep,*,end=111,err=111)
470 & (((ffthet(llll,lll,ll,i,j,k,iblock),
471 & ffthet(lll,llll,ll,i,j,k,iblock),
472 & ggthet(llll,lll,ll,i,j,k,iblock)
473 & ,ggthet(lll,llll,ll,i,j,k,iblock),
474 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
478 C write(iout,*) "KURWA1.1"
480 C For dummy ends assign glycine-type coefficients of theta-only terms; the
481 C coefficients of theta-and-gamma-dependent terms are zero.
486 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
487 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
489 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
490 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
493 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
495 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
498 C write(iout,*) "KURWA1.5"
499 C Substitution for D aminoacids from symmetry.
502 do j=-nthetyp,nthetyp
503 do k=-nthetyp,nthetyp
504 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
506 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
510 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
511 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
512 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
513 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
519 ffthet(llll,lll,ll,i,j,k,iblock)=
520 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
521 ffthet(lll,llll,ll,i,j,k,iblock)=
522 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
523 ggthet(llll,lll,ll,i,j,k,iblock)=
524 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
525 ggthet(lll,llll,ll,i,j,k,iblock)=
526 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
536 C Control printout of the coefficients of virtual-bond-angle potentials
539 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
544 write (iout,'(//4a)')
545 & 'Type ',onelett(i),onelett(j),onelett(k)
546 write (iout,'(//a,10x,a)') " l","a[l]"
547 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
548 write (iout,'(i2,1pe15.5)')
549 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
551 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
552 & "b",l,"c",l,"d",l,"e",l
554 write (iout,'(i2,4(1pe15.5))') m,
555 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
556 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
560 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
561 & "f+",l,"f-",l,"g+",l,"g-",l
564 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
565 & ffthet(n,m,l,i,j,k,iblock),
566 & ffthet(m,n,l,i,j,k,iblock),
567 & ggthet(n,m,l,i,j,k,iblock),
568 & ggthet(m,n,l,i,j,k,iblock)
581 C here will be the apropriate recalibrating for D-aminoacid
582 read (ithep,*,end=111,err=111) nthetyp
583 do i=-nthetyp+1,nthetyp-1
584 read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
585 do j=0,nbend_kcc_Tb(i)
586 read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
591 & "Parameters of the valence-only potentials"
592 do i=-nthetyp+1,nthetyp-1
593 write (iout,'(2a)') "Type ",toronelet(i)
594 do k=0,nbend_kcc_Tb(i)
595 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
604 C write(iout,*) 'KURWA2'
607 C Read the parameters of the probability distribution/energy expression
608 C of the side chains.
611 cc write (iout,*) "tu dochodze",i
612 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
616 dsc_inv(i)=1.0D0/dsc(i)
627 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
628 & ((blower(k,l,1),l=1,k),k=1,3)
629 censc(1,1,-i)=censc(1,1,i)
630 censc(2,1,-i)=censc(2,1,i)
631 censc(3,1,-i)=-censc(3,1,i)
633 read (irotam,*,end=112,err=112) bsc(j,i)
634 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
635 & ((blower(k,l,j),l=1,k),k=1,3)
636 censc(1,j,-i)=censc(1,j,i)
637 censc(2,j,-i)=censc(2,j,i)
638 censc(3,j,-i)=-censc(3,j,i)
639 C BSC is amplitude of Gaussian
646 akl=akl+blower(k,m,j)*blower(l,m,j)
650 if (((k.eq.3).and.(l.ne.3))
651 & .or.((l.eq.3).and.(k.ne.3))) then
652 gaussc(k,l,j,-i)=-akl
653 gaussc(l,k,j,-i)=-akl
665 write (iout,'(/a)') 'Parameters of side-chain local geometry'
669 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
670 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
671 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
672 c write (iout,'(a,f10.4,4(16x,f10.4))')
673 c & 'Center ',(bsc(j,i),j=1,nlobi)
674 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
675 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
676 & 'log h',(bsc(j,i),j=1,nlobi)
677 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
678 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
685 c blower(k,l,j)=gaussc(ind,j,i)
690 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
691 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
698 C Read scrot parameters for potentials determined from all-atom AM1 calculations
699 C added by Urszula Kozlowska 07/11/2007
702 read (irotam,*,end=112,err=112)
704 read (irotam,*,end=112,err=112)
707 read(irotam,*,end=112,err=112) sc_parmin(j,i)
714 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
715 C interaction energy of the Gly, Ala, and Pro prototypes.
717 read (ifourier,*,end=115,err=115) nloctyp
718 SPLIT_FOURIERTOR = nloctyp.lt.0
719 nloctyp = iabs(nloctyp)
721 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
722 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
723 itype2loc(ntyp1)=nloctyp
724 iloctyp(nloctyp)=ntyp1
726 itype2loc(-i)=-itype2loc(i)
735 iloctyp(-i)=-iloctyp(i)
737 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
738 c write (iout,*) "nloctyp",nloctyp,
739 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
742 c write (iout,*) "NEWCORR",i
743 read (ifourier,*,end=115,err=115)
746 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
749 c write (iout,*) "NEWCORR BNEW1"
750 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
753 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
756 c write (iout,*) "NEWCORR BNEW2"
757 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
759 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
760 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
762 c write (iout,*) "NEWCORR CCNEW"
763 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
765 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
766 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
768 c write (iout,*) "NEWCORR DDNEW"
769 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
773 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
777 c write (iout,*) "NEWCORR EENEW1"
778 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
780 read (ifourier,*,end=115,err=115) e0new(ii,i)
782 c write (iout,*) (e0new(ii,i),ii=1,3)
784 c write (iout,*) "NEWCORR EENEW"
787 ccnew(ii,1,i)=ccnew(ii,1,i)/2
788 ccnew(ii,2,i)=ccnew(ii,2,i)/2
789 ddnew(ii,1,i)=ddnew(ii,1,i)/2
790 ddnew(ii,2,i)=ddnew(ii,2,i)/2
795 bnew1(ii,1,-i)= bnew1(ii,1,i)
796 bnew1(ii,2,-i)=-bnew1(ii,2,i)
797 bnew2(ii,1,-i)= bnew2(ii,1,i)
798 bnew2(ii,2,-i)=-bnew2(ii,2,i)
801 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
802 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
803 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
804 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
805 ccnew(ii,1,-i)=ccnew(ii,1,i)
806 ccnew(ii,2,-i)=-ccnew(ii,2,i)
807 ddnew(ii,1,-i)=ddnew(ii,1,i)
808 ddnew(ii,2,-i)=-ddnew(ii,2,i)
810 e0new(1,-i)= e0new(1,i)
811 e0new(2,-i)=-e0new(2,i)
812 e0new(3,-i)=-e0new(3,i)
814 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
815 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
816 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
817 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
821 write (iout,'(a)') "Coefficients of the multibody terms"
822 do i=-nloctyp+1,nloctyp-1
823 write (iout,*) "Type: ",onelet(iloctyp(i))
824 write (iout,*) "Coefficients of the expansion of B1"
826 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
828 write (iout,*) "Coefficients of the expansion of B2"
830 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
832 write (iout,*) "Coefficients of the expansion of C"
833 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
834 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
835 write (iout,*) "Coefficients of the expansion of D"
836 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
837 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
838 write (iout,*) "Coefficients of the expansion of E"
839 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
842 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
847 IF (SPLIT_FOURIERTOR) THEN
849 c write (iout,*) "NEWCORR TOR",i
850 read (ifourier,*,end=115,err=115)
853 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
856 c write (iout,*) "NEWCORR BNEW1 TOR"
857 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
860 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
863 c write (iout,*) "NEWCORR BNEW2 TOR"
864 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
866 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
867 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
869 c write (iout,*) "NEWCORR CCNEW TOR"
870 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
872 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
873 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
875 c write (iout,*) "NEWCORR DDNEW TOR"
876 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
880 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
884 c write (iout,*) "NEWCORR EENEW1 TOR"
885 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
887 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
889 c write (iout,*) (e0newtor(ii,i),ii=1,3)
891 c write (iout,*) "NEWCORR EENEW TOR"
894 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
895 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
896 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
897 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
902 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
903 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
904 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
905 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
908 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
909 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
910 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
911 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
913 e0newtor(1,-i)= e0newtor(1,i)
914 e0newtor(2,-i)=-e0newtor(2,i)
915 e0newtor(3,-i)=-e0newtor(3,i)
917 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
918 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
919 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
920 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
925 & "Single-body coefficients of the torsional potentials"
926 do i=-nloctyp+1,nloctyp-1
927 write (iout,*) "Type: ",onelet(iloctyp(i))
928 write (iout,*) "Coefficients of the expansion of B1tor"
930 write (iout,'(3hB1(,i1,1h),3f10.5)')
931 & j,(bnew1tor(k,j,i),k=1,3)
933 write (iout,*) "Coefficients of the expansion of B2tor"
935 write (iout,'(3hB2(,i1,1h),3f10.5)')
936 & j,(bnew2tor(k,j,i),k=1,3)
938 write (iout,*) "Coefficients of the expansion of Ctor"
939 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
940 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
941 write (iout,*) "Coefficients of the expansion of Dtor"
942 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
943 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
944 write (iout,*) "Coefficients of the expansion of Etor"
945 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
948 write (iout,'(1hE,2i1,2f10.5)')
949 & j,k,(eenewtor(l,j,k,i),l=1,2)
955 do i=-nloctyp+1,nloctyp-1
958 bnew1tor(ii,j,i)=bnew1(ii,j,i)
963 bnew2tor(ii,j,i)=bnew2(ii,j,i)
967 ccnewtor(ii,1,i)=ccnew(ii,1,i)
968 ccnewtor(ii,2,i)=ccnew(ii,2,i)
969 ddnewtor(ii,1,i)=ddnew(ii,1,i)
970 ddnewtor(ii,2,i)=ddnew(ii,2,i)
973 ENDIF !SPLIT_FOURIER_TOR
976 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
978 read (ifourier,*,end=115,err=115)
979 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
981 write (iout,*) 'Type ',onelet(iloctyp(i))
982 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
996 c B1tilde(1,i) = b(3)
997 c B1tilde(2,i) =-b(5)
998 c B1tilde(1,-i) =-b(3)
999 c B1tilde(2,-i) =b(5)
1000 c b1tilde(1,i)=0.0d0
1001 c b1tilde(2,i)=0.0d0
1006 cc B1tilde(1,i) = b(3,i)
1007 cc B1tilde(2,i) =-b(5,i)
1008 C B1tilde(1,-i) =-b(3,i)
1009 C B1tilde(2,-i) =b(5,i)
1010 cc b1tilde(1,i)=0.0d0
1011 cc b1tilde(2,i)=0.0d0
1019 CCold(1,1,i)= b(7,i)
1020 CCold(2,2,i)=-b(7,i)
1021 CCold(2,1,i)= b(9,i)
1022 CCold(1,2,i)= b(9,i)
1023 CCold(1,1,-i)= b(7,i)
1024 CCold(2,2,-i)=-b(7,i)
1025 CCold(2,1,-i)=-b(9,i)
1026 CCold(1,2,-i)=-b(9,i)
1031 c Ctilde(1,1,i)= CCold(1,1,i)
1032 c Ctilde(1,2,i)= CCold(1,2,i)
1033 c Ctilde(2,1,i)=-CCold(2,1,i)
1034 c Ctilde(2,2,i)=-CCold(2,2,i)
1036 c Ctilde(1,1,i)=0.0d0
1037 c Ctilde(1,2,i)=0.0d0
1038 c Ctilde(2,1,i)=0.0d0
1039 c Ctilde(2,2,i)=0.0d0
1040 DDold(1,1,i)= b(6,i)
1041 DDold(2,2,i)=-b(6,i)
1042 DDold(2,1,i)= b(8,i)
1043 DDold(1,2,i)= b(8,i)
1044 DDold(1,1,-i)= b(6,i)
1045 DDold(2,2,-i)=-b(6,i)
1046 DDold(2,1,-i)=-b(8,i)
1047 DDold(1,2,-i)=-b(8,i)
1052 c Dtilde(1,1,i)= DD(1,1,i)
1053 c Dtilde(1,2,i)= DD(1,2,i)
1054 c Dtilde(2,1,i)=-DD(2,1,i)
1055 c Dtilde(2,2,i)=-DD(2,2,i)
1057 c Dtilde(1,1,i)=0.0d0
1058 c Dtilde(1,2,i)=0.0d0
1059 c Dtilde(2,1,i)=0.0d0
1060 c Dtilde(2,2,i)=0.0d0
1061 EEold(1,1,i)= b(10,i)+b(11,i)
1062 EEold(2,2,i)=-b(10,i)+b(11,i)
1063 EEold(2,1,i)= b(12,i)-b(13,i)
1064 EEold(1,2,i)= b(12,i)+b(13,i)
1065 EEold(1,1,-i)= b(10,i)+b(11,i)
1066 EEold(2,2,-i)=-b(10,i)+b(11,i)
1067 EEold(2,1,-i)=-b(12,i)+b(13,i)
1068 EEold(1,2,-i)=-b(12,i)-b(13,i)
1069 write(iout,*) "TU DOCHODZE"
1075 c ee(2,1,i)=ee(1,2,i)
1080 &"Coefficients of the cumulants (independent of valence angles)"
1081 do i=-nloctyp+1,nloctyp-1
1082 write (iout,*) 'Type ',onelet(iloctyp(i))
1084 write(iout,'(2f10.5)') B(3,i),B(5,i)
1086 write(iout,'(2f10.5)') B(2,i),B(4,i)
1089 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1093 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1097 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1102 C write (iout,*) 'KURWAKURWA'
1105 C Read torsional parameters in old format
1107 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1108 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1109 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1112 read (itorp,'(a)',end=113,err=113)
1114 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1120 write (iout,'(/a/)') 'Torsional constants:'
1123 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1124 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1130 C Read torsional parameters
1132 IF (TOR_MODE.eq.0) THEN
1134 read (itorp,*,end=113,err=113) ntortyp
1135 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1137 itype2loc(i)=itortyp(i)
1139 write (iout,*) 'ntortyp',ntortyp
1141 itype2loc(-i)=-itype2loc(i)
1143 itortyp(ntyp1)=ntortyp
1146 itortyp(i)=-itortyp(-i)
1148 c write (iout,*) 'ntortyp',ntortyp
1150 do j=-ntortyp+1,ntortyp-1
1151 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1153 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1154 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1157 do k=1,nterm(i,j,iblock)
1158 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1160 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1161 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1162 v0ij=v0ij+si*v1(k,i,j,iblock)
1165 do k=1,nlor(i,j,iblock)
1166 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1167 & vlor2(k,i,j),vlor3(k,i,j)
1168 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1171 v0(-i,-j,iblock)=v0ij
1177 write (iout,'(/a/)') 'Torsional constants:'
1180 write (iout,*) 'ityp',i,' jtyp',j
1181 write (iout,*) 'Fourier constants'
1182 do k=1,nterm(i,j,iblock)
1183 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1186 write (iout,*) 'Lorenz constants'
1187 do k=1,nlor(i,j,iblock)
1188 write (iout,'(3(1pe15.5))')
1189 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1195 C 6/23/01 Read parameters for double torsionals
1199 do j=-ntortyp+1,ntortyp-1
1200 do k=-ntortyp+1,ntortyp-1
1201 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1202 c write (iout,*) "OK onelett",
1205 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1206 & .or. t3.ne.toronelet(k)) then
1207 write (iout,*) "Error in double torsional parameter file",
1210 call MPI_Finalize(Ierror)
1212 stop "Error in double torsional parameter file"
1214 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1215 & ntermd_2(i,j,k,iblock)
1216 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1217 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1218 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1219 & ntermd_1(i,j,k,iblock))
1220 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1221 & ntermd_1(i,j,k,iblock))
1222 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1223 & ntermd_1(i,j,k,iblock))
1224 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1225 & ntermd_1(i,j,k,iblock))
1226 C Martix of D parameters for one dimesional foureir series
1227 do l=1,ntermd_1(i,j,k,iblock)
1228 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1229 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1230 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1231 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1232 c write(iout,*) "whcodze" ,
1233 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1235 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1236 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1237 & v2s(m,l,i,j,k,iblock),
1238 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1239 C Martix of D parameters for two dimesional fourier series
1240 do l=1,ntermd_2(i,j,k,iblock)
1242 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1243 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1244 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1245 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1254 write (iout,*) 'Constants for double torsionals'
1257 do j=-ntortyp+1,ntortyp-1
1258 do k=-ntortyp+1,ntortyp-1
1259 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1260 & ' nsingle',ntermd_1(i,j,k,iblock),
1261 & ' ndouble',ntermd_2(i,j,k,iblock)
1263 write (iout,*) 'Single angles:'
1264 do l=1,ntermd_1(i,j,k,iblock)
1265 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1266 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1267 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1268 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1271 write (iout,*) 'Pairs of angles:'
1272 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1273 do l=1,ntermd_2(i,j,k,iblock)
1274 write (iout,'(i5,20f10.5)')
1275 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1278 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1279 do l=1,ntermd_2(i,j,k,iblock)
1280 write (iout,'(i5,20f10.5)')
1281 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1282 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1291 ELSE IF (TOR_MODE.eq.1) THEN
1293 C read valence-torsional parameters
1294 read (itorp,*,end=113,err=113) ntortyp
1296 write (iout,*) "Valence-torsional parameters read in ntortyp",
1298 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1299 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1301 itortyp(i)=-itortyp(-i)
1303 do i=-ntortyp+1,ntortyp-1
1304 do j=-ntortyp+1,ntortyp-1
1305 C first we read the cos and sin gamma parameters
1306 read (itorp,'(13x,a)',end=113,err=113) string
1307 write (iout,*) i,j,string
1308 read (itorp,*,end=113,err=113)
1309 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1310 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1311 do k=1,nterm_kcc(j,i)
1312 do l=1,nterm_kcc_Tb(j,i)
1313 do ll=1,nterm_kcc_Tb(j,i)
1314 read (itorp,*,end=113,err=113) ii,jj,kk,
1315 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1323 c AL 4/8/16: Calculate coefficients from one-body parameters
1326 itortyp(i)=itype2loc(i)
1328 if (lprint) write (iout,*)
1329 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1331 do i=-ntortyp+1,ntortyp-1
1332 do j=-ntortyp+1,ntortyp-1
1335 do k=1,nterm_kcc_Tb(j,i)
1336 do l=1,nterm_kcc_Tb(j,i)
1337 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1338 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1339 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1340 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1343 do k=1,nterm_kcc_Tb(j,i)
1344 do l=1,nterm_kcc_Tb(j,i)
1346 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1347 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1348 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1349 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1351 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1352 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1353 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1354 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1358 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)
1364 if (tor_mode.gt.0 .and. lprint) then
1365 c Print valence-torsional parameters
1367 & "Parameters of the valence-torsional potentials"
1368 do i=-ntortyp+1,ntortyp-1
1369 do j=-ntortyp+1,ntortyp-1
1370 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1371 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1372 do k=1,nterm_kcc(j,i)
1373 do l=1,nterm_kcc_Tb(j,i)
1374 do ll=1,nterm_kcc_Tb(j,i)
1375 write (iout,'(3i5,2f15.4)')
1376 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1385 C Read of Side-chain backbone correlation parameters
1386 C Modified 11 May 2012 by Adasko
1389 read (isccor,*,end=119,err=119) nsccortyp
1390 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1392 isccortyp(i)=-isccortyp(-i)
1394 iscprol=isccortyp(20)
1395 c write (iout,*) 'ntortyp',ntortyp
1397 cc maxinter is maximum interaction sites
1401 read (isccor,*,end=119,err=119)
1402 &nterm_sccor(i,j),nlor_sccor(i,j)
1403 c write (iout,*) nterm_sccor(i,j)
1409 nterm_sccor(-i,j)=nterm_sccor(i,j)
1410 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1411 nterm_sccor(i,-j)=nterm_sccor(i,j)
1412 c write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1413 c & nterm_sccor(-i,-j),nterm_sccor(i,-j)
1414 do k=1,nterm_sccor(i,j)
1415 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1417 if (j.eq.iscprol) then
1418 if (i.eq.isccortyp(10)) then
1419 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1420 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1422 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1423 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1424 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1425 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1426 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1427 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1428 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1429 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1432 if (i.eq.isccortyp(10)) then
1433 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1434 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1436 if (j.eq.isccortyp(10)) then
1437 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1438 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1440 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1441 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1442 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1443 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1444 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1445 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1449 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1450 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1451 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1452 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1455 do k=1,nlor_sccor(i,j)
1456 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1457 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1458 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1459 &(1+vlor3sccor(k,i,j)**2)
1461 v0sccor(l,i,j)=v0ijsccor
1462 v0sccor(l,-i,j)=v0ijsccor1
1463 v0sccor(l,i,-j)=v0ijsccor2
1464 v0sccor(l,-i,-j)=v0ijsccor3
1470 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1474 write (iout,*) 'ityp',i,' jtyp',j
1475 write (iout,*) 'Fourier constants'
1476 do k=1,nterm_sccor(i,j)
1477 write (iout,'(2(1pe15.5))')
1478 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1480 write (iout,*) 'Lorenz constants'
1481 do k=1,nlor_sccor(i,j)
1482 write (iout,'(3(1pe15.5))')
1483 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1490 C Read electrostatic-interaction parameters
1493 write (iout,'(/a)') 'Electrostatic interaction constants:'
1494 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1495 & 'IT','JT','APP','BPP','AEL6','AEL3'
1497 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1498 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1499 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1500 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1505 app (i,j)=epp(i,j)*rri*rri
1506 bpp (i,j)=-2.0D0*epp(i,j)*rri
1507 ael6(i,j)=elpp6(i,j)*4.2D0**6
1508 ael3(i,j)=elpp3(i,j)*4.2D0**3
1509 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1510 & ael6(i,j),ael3(i,j)
1514 C Read side-chain interaction parameters.
1516 read (isidep,*,end=117,err=117) ipot,expon
1517 if (ipot.lt.1 .or. ipot.gt.5) then
1518 write (iout,'(2a)') 'Error while reading SC interaction',
1519 & 'potential file - unknown potential type.'
1523 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1524 & ', exponents are ',expon,2*expon
1525 goto (10,20,30,30,40) ipot
1526 C----------------------- LJ potential ---------------------------------
1527 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1528 & (sigma0(i),i=1,ntyp)
1530 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1531 write (iout,'(a/)') 'The epsilon array:'
1532 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1533 write (iout,'(/a)') 'One-body parameters:'
1534 write (iout,'(a,4x,a)') 'residue','sigma'
1535 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1538 C----------------------- LJK potential --------------------------------
1539 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1540 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1542 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1543 write (iout,'(a/)') 'The epsilon array:'
1544 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1545 write (iout,'(/a)') 'One-body parameters:'
1546 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1547 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1551 C---------------------- GB or BP potential -----------------------------
1553 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1555 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1556 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1557 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1558 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1560 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1561 C write(iout,*) "WARNING!!",i,ntyp
1562 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1564 C epslip(i,j)=epslip(i,j)+0.05d0
1567 C For the GB potential convert sigma'**2 into chi'
1570 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1574 write (iout,'(/a/)') 'Parameters of the BP potential:'
1575 write (iout,'(a/)') 'The epsilon array:'
1576 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1577 write (iout,'(/a)') 'One-body parameters:'
1578 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1580 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1581 & chip(i),alp(i),i=1,ntyp)
1584 C--------------------- GBV potential -----------------------------------
1585 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1586 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1587 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1589 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1590 write (iout,'(a/)') 'The epsilon array:'
1591 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1592 write (iout,'(/a)') 'One-body parameters:'
1593 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1594 & 's||/s_|_^2',' chip ',' alph '
1595 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1596 & sigii(i),chip(i),alp(i),i=1,ntyp)
1600 C-----------------------------------------------------------------------
1601 C Calculate the "working" parameters of SC interactions.
1605 epslip(i,j)=epslip(j,i)
1610 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1611 sigma(j,i)=sigma(i,j)
1612 rs0(i,j)=dwa16*sigma(i,j)
1616 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1617 & 'Working parameters of the SC interactions:',
1618 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1623 epsijlip=epslip(i,j)
1624 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1633 sigeps=dsign(1.0D0,epsij)
1635 aa_aq(i,j)=epsij*rrij*rrij
1636 bb_aq(i,j)=-sigeps*epsij*rrij
1637 aa_aq(j,i)=aa_aq(i,j)
1638 bb_aq(j,i)=bb_aq(i,j)
1639 sigeps=dsign(1.0D0,epsijlip)
1640 epsijlip=dabs(epsijlip)
1641 aa_lip(i,j)=epsijlip*rrij*rrij
1642 bb_lip(i,j)=-sigeps*epsijlip*rrij
1643 aa_lip(j,i)=aa_lip(i,j)
1644 bb_lip(j,i)=bb_lip(i,j)
1646 sigt1sq=sigma0(i)**2
1647 sigt2sq=sigma0(j)**2
1650 ratsig1=sigt2sq/sigt1sq
1651 ratsig2=1.0D0/ratsig1
1652 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1653 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1654 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1658 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1659 sigmaii(i,j)=rsum_max
1660 sigmaii(j,i)=rsum_max
1662 c sigmaii(i,j)=r0(i,j)
1663 c sigmaii(j,i)=r0(i,j)
1665 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1666 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1667 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1668 augm(i,j)=epsij*r_augm**(2*expon)
1669 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1676 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1677 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1678 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1683 C Define the SC-p interaction constants
1687 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1689 c aad(i,1)=0.3D0*4.0D0**12
1690 C Following line for constants currently implemented
1691 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1692 aad(i,1)=1.5D0*4.0D0**12
1693 c aad(i,1)=0.17D0*5.6D0**12
1695 C "Soft" SC-p repulsion
1697 C Following line for constants currently implemented
1698 c aad(i,1)=0.3D0*4.0D0**6
1699 C "Hard" SC-p repulsion
1700 bad(i,1)=3.0D0*4.0D0**6
1701 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1710 C 8/9/01 Read the SC-p interaction constants from file
1713 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1716 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1717 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1718 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1719 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1723 write (iout,'(/a)') "Parameters of SC-p interactions:"
1725 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1726 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1731 C Define the constants of the disulfide bridge
1735 c Old arbitrary potential - commented out.
1740 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1741 c energy surface of diethyl disulfide.
1742 c A. Liwo and U. Kozlowska, 11/24/03
1752 & write (iout,*) 'Dynamic formation/breaking of disulfides'
1754 ss_depth=ebr/wsc-0.25*eps(1,1)
1755 C write(iout,*) akcm,whpb,wsc,'KURWA'
1756 Ht=Ht/wsc-0.25*eps(1,1)
1765 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1770 write (iout,'(/a)') "Disulfide bridge parameters:"
1771 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1772 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1773 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1774 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1776 write (iout,'(a)') "Parameters of the 'trisulfide' potential"
1777 write(iout,*) "ATRISS",atriss
1778 write(iout,*) "BTRISS",btriss
1779 write(iout,*) "CTRISS",ctriss
1780 write(iout,*) "DTRISS",dtriss
1783 if (shield_mode.gt.0) then
1785 C VSolvSphere the volume of solving sphere
1787 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1788 C there will be no distinction between proline peptide group and normal peptide
1789 C group in case of shielding parameters
1790 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1791 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1792 write (iout,*) "VSolvSphere",VSolvSphere," VSolvSphere_div",
1794 C long axis of side chain
1796 long_r_sidechain(i)=vbldsc0(1,i)
1797 short_r_sidechain(i)=sigma0(i)
1802 111 write (iout,*) "Error reading bending energy parameters."
1804 112 write (iout,*) "Error reading rotamer energy parameters."
1806 113 write (iout,*) "Error reading torsional energy parameters."
1808 114 write (iout,*) "Error reading double torsional energy parameters."
1811 & "Error reading cumulant (multibody energy) parameters."
1813 116 write (iout,*) "Error reading electrostatic energy parameters."
1815 1161 write (iout,*) "Error reading lipid parameters."
1817 117 write (iout,*) "Error reading side chain interaction parameters."
1819 118 write (iout,*) "Error reading SCp interaction parameters."
1821 119 write (iout,*) "Error reading SCCOR parameters"
1823 121 write (iout,*) "Error reading bond parameters"
1826 call MPI_Finalize(Ierror)