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 r0_corr=cutoff_corr-delt_corr
68 write (iout,*) "iparm",iparm," myparm",myparm
69 c If reading not own parameters, skip assignment
70 call reada(controlcard,"D0CM",d0cm,3.78d0)
71 call reada(controlcard,"AKCM",akcm,15.1d0)
72 call reada(controlcard,"AKTH",akth,11.0d0)
73 call reada(controlcard,"AKCT",akct,12.0d0)
74 call reada(controlcard,"V1SS",v1ss,-1.08d0)
75 call reada(controlcard,"V2SS",v2ss,7.61d0)
76 call reada(controlcard,"V3SS",v3ss,13.7d0)
77 call reada(controlcard,"EBR",ebr,-5.50D0)
78 call reada(controlcard,"DTRISS",dtriss,1.0D0)
79 call reada(controlcard,"ATRISS",atriss,0.3D0)
80 call reada(controlcard,"BTRISS",btriss,0.02D0)
81 call reada(controlcard,"CTRISS",ctriss,1.0D0)
82 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
84 C dyn_ss_mask(i)=.false.
88 c Old arbitrary potential - commented out.
93 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
94 c energy surface of diethyl disulfide.
95 c A. Liwo and U. Kozlowska, 11/24/03
107 dyn_ssbond_ij(i,j)=1.0d300
110 call reada(controlcard,"HT",Ht,0.0D0)
112 C ss_depth=ebr/wsc-0.25*eps(1,1)
113 C write(iout,*) HT,wsc,eps(1,1),'KURWA'
114 C Ht=Ht/wsc-0.25*eps(1,1)
123 C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
126 if (iparm.eq.myparm .or. .not.separate_parset) then
129 c Setup weights for UNRES
153 c write (iout,*) "PARMREAD: wsaxs",wsaxs
158 c write(iout,*)"PARMREAD: wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
159 c & " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
163 C Juyong:READ init_vars
164 C Initialize variables!
165 C Juyong:READ read_info
166 C READ fragment information!!
167 C both routines should be in dfa.F file!!
168 write (iout,*) "Before initializing DFA",wdfa_dist,wdfa_tor,
170 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
171 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
172 write (iout,*) "Calling init_dfa_vars"
175 write (iout,*) 'init_dfa_vars finished!'
178 write (iout,*) 'read_dfa_info finished!'
182 call card_concat(controlcard,.false.)
184 c Return if not own parameters
186 if (iparm.ne.myparm .and. separate_parset) return
188 call reads(controlcard,"BONDPAR",bondname_t,bondname)
189 open (ibond,file=bondname_t,status='old')
191 call reads(controlcard,"THETPAR",thetname_t,thetname)
192 open (ithep,file=thetname_t,status='old')
194 call reads(controlcard,"ROTPAR",rotname_t,rotname)
195 open (irotam,file=rotname_t,status='old')
197 call reads(controlcard,"TORPAR",torname_t,torname)
198 open (itorp,file=torname_t,status='old')
200 call reads(controlcard,"TORDPAR",tordname_t,tordname)
201 write (iout,*) "tor_mode",tor_mode
204 & open (itordp,file=tordname_t,status='old')
206 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
207 open (isccor,file=sccorname_t,status='old')
209 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
210 open (ifourier,file=fouriername_t,status='old')
212 call reads(controlcard,"ELEPAR",elename_t,elename)
213 open (ielep,file=elename_t,status='old')
215 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
216 open (isidep,file=sidename_t,status='old')
218 call reads(controlcard,"SCPPAR",scpname_t,scpname)
219 open (iscpp,file=scpname_t,status='old')
221 write (iout,*) "Parameter set:",iparm
222 write (iout,*) "Energy-term weights:"
224 write (iout,'(a16,f10.5)') wname(i),ww(i)
226 write (iout,*) "Sidechain potential file : ",
227 & sidename_t(:ilen(sidename_t))
229 write (iout,*) "SCp potential file : ",
230 & scpname_t(:ilen(scpname_t))
232 write (iout,*) "Electrostatic potential file : ",
233 & elename_t(:ilen(elename_t))
234 write (iout,*) "Cumulant coefficient file : ",
235 & fouriername_t(:ilen(fouriername_t))
236 write (iout,*) "Torsional parameter file : ",
237 & torname_t(:ilen(torname_t))
238 write (iout,*) "Double torsional parameter file : ",
239 & tordname_t(:ilen(tordname_t))
240 write (iout,*) "Backbone-rotamer parameter file : ",
241 & sccorname(:ilen(sccorname))
242 write (iout,*) "Bond & inertia constant file : ",
243 & bondname_t(:ilen(bondname_t))
244 write (iout,*) "Bending parameter file : ",
245 & thetname_t(:ilen(thetname_t))
246 write (iout,*) "Rotamer parameter file : ",
247 & rotname_t(:ilen(rotname_t))
250 c Read the virtual-bond parameters, masses, and moments of inertia
251 c and Stokes' radii of the peptide group and side chains
254 read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp,mp,ip,pstok
257 read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i),
258 & msc(i),isc(i),restok(i)
259 dsc(i) = vbldsc0(1,i)
263 dsc_inv(i)=1.0D0/dsc(i)
267 read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk,
270 read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
271 & aksc(j,i),abond0(j,i),j=1,nbondterm(i)),msc(i),isc(i),restok(i)
272 dsc(i) = vbldsc0(1,i)
276 dsc_inv(i)=1.0D0/dsc(i)
281 write(iout,'(/a/)')"Force constants virtual bonds:"
282 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
284 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
286 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
287 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
289 write (iout,'(13x,3f10.5)')
290 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
294 c write (iout,*) "iliptranpar",iliptranpar
295 c write (iout,*) "liptranname ",liptranname
296 read(iliptranpar,*,end=1161,err=1161) pepliptran
298 read(iliptranpar,*,end=1161,err=1161) liptranene(i)
302 write (iout,'(/a)') "Water-lipid transfer parameters"
303 write (iout,'(a3,3x,f10.5)') 'p',pepliptran
305 write (iout,'(a3,3x,f10.5)') restyp(i),liptranene(i)
310 C Read the parameters of the probability distribution/energy expression
311 C of the virtual-bond valence angles theta
314 read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
315 & (bthet(j,i,1,1),j=1,2)
316 read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
317 read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
318 read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
322 athet(1,i,1,-1)=athet(1,i,1,1)
323 athet(2,i,1,-1)=athet(2,i,1,1)
324 bthet(1,i,1,-1)=-bthet(1,i,1,1)
325 bthet(2,i,1,-1)=-bthet(2,i,1,1)
326 athet(1,i,-1,1)=-athet(1,i,1,1)
327 athet(2,i,-1,1)=-athet(2,i,1,1)
328 bthet(1,i,-1,1)=bthet(1,i,1,1)
329 bthet(2,i,-1,1)=bthet(2,i,1,1)
333 athet(1,i,-1,-1)=athet(1,-i,1,1)
334 athet(2,i,-1,-1)=-athet(2,-i,1,1)
335 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
336 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
337 athet(1,i,-1,1)=athet(1,-i,1,1)
338 athet(2,i,-1,1)=-athet(2,-i,1,1)
339 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
340 bthet(2,i,-1,1)=bthet(2,-i,1,1)
341 athet(1,i,1,-1)=-athet(1,-i,1,1)
342 athet(2,i,1,-1)=athet(2,-i,1,1)
343 bthet(1,i,1,-1)=bthet(1,-i,1,1)
344 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
349 polthet(j,i)=polthet(j,-i)
352 gthet(j,i)=gthet(j,-i)
358 c & 'Parameters of the virtual-bond valence angles:'
359 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
360 c & ' ATHETA0 ',' A1 ',' A2 ',
363 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
364 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
366 c write (iout,'(/a/9x,5a/79(1h-))')
367 c & 'Parameters of the expression for sigma(theta_c):',
368 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
369 c & ' ALPH3 ',' SIGMA0C '
371 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
372 c & (polthet(j,i),j=0,3),sigc0(i)
374 c write (iout,'(/a/9x,5a/79(1h-))')
375 c & 'Parameters of the second gaussian:',
376 c & ' THETA0 ',' SIGMA0 ',' G1 ',
379 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
380 c & sig0(i),(gthet(j,i),j=1,3)
383 & 'Parameters of the virtual-bond valence angles:'
384 write (iout,'(/a/9x,5a/79(1h-))')
385 & 'Coefficients of expansion',
386 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
387 & ' b1*10^1 ',' b2*10^1 '
389 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
390 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
391 & (10*bthet(j,i,1,1),j=1,2)
393 write (iout,'(/a/9x,5a/79(1h-))')
394 & 'Parameters of the expression for sigma(theta_c):',
395 & ' alpha0 ',' alph1 ',' alph2 ',
396 & ' alhp3 ',' sigma0c '
398 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
399 & (polthet(j,i),j=0,3),sigc0(i)
401 write (iout,'(/a/9x,5a/79(1h-))')
402 & 'Parameters of the second gaussian:',
403 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
406 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
407 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
411 IF (tor_mode.eq.0) THEN
413 C Read the parameters of Utheta determined from ab initio surfaces
414 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
416 read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
417 & ntheterm3,nsingle,ndouble
418 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
419 read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
421 ithetyp(i)=-ithetyp(-i)
424 do i=-maxthetyp,maxthetyp
425 do j=-maxthetyp,maxthetyp
426 do k=-maxthetyp,maxthetyp
427 aa0thet(i,j,k,iblock)=0.0d0
429 aathet(l,i,j,k,iblock)=0.0d0
433 bbthet(m,l,i,j,k,iblock)=0.0d0
434 ccthet(m,l,i,j,k,iblock)=0.0d0
435 ddthet(m,l,i,j,k,iblock)=0.0d0
436 eethet(m,l,i,j,k,iblock)=0.0d0
442 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
443 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
451 C write (iout,*) "KURWA1"
454 do j=-nthetyp,nthetyp
455 do k=-nthetyp,nthetyp
456 read (ithep,'(6a)',end=111,err=111) res1
457 write(iout,*) res1,i,j,k
458 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
459 read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
461 read (ithep,*,end=111,err=111)
462 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
463 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
464 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
465 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
467 read (ithep,*,end=111,err=111)
468 & (((ffthet(llll,lll,ll,i,j,k,iblock),
469 & ffthet(lll,llll,ll,i,j,k,iblock),
470 & ggthet(llll,lll,ll,i,j,k,iblock)
471 & ,ggthet(lll,llll,ll,i,j,k,iblock),
472 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
476 C write(iout,*) "KURWA1.1"
478 C For dummy ends assign glycine-type coefficients of theta-only terms; the
479 C coefficients of theta-and-gamma-dependent terms are zero.
484 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
485 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
487 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
488 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
491 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
493 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
496 C write(iout,*) "KURWA1.5"
497 C Substitution for D aminoacids from symmetry.
500 do j=-nthetyp,nthetyp
501 do k=-nthetyp,nthetyp
502 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
504 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
508 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
509 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
510 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
511 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
517 ffthet(llll,lll,ll,i,j,k,iblock)=
518 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
519 ffthet(lll,llll,ll,i,j,k,iblock)=
520 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
521 ggthet(llll,lll,ll,i,j,k,iblock)=
522 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
523 ggthet(lll,llll,ll,i,j,k,iblock)=
524 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
534 C Control printout of the coefficients of virtual-bond-angle potentials
537 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
542 write (iout,'(//4a)')
543 & 'Type ',onelett(i),onelett(j),onelett(k)
544 write (iout,'(//a,10x,a)') " l","a[l]"
545 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
546 write (iout,'(i2,1pe15.5)')
547 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
549 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
550 & "b",l,"c",l,"d",l,"e",l
552 write (iout,'(i2,4(1pe15.5))') m,
553 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
554 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
558 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
559 & "f+",l,"f-",l,"g+",l,"g-",l
562 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
563 & ffthet(n,m,l,i,j,k,iblock),
564 & ffthet(m,n,l,i,j,k,iblock),
565 & ggthet(n,m,l,i,j,k,iblock),
566 & ggthet(m,n,l,i,j,k,iblock)
579 C here will be the apropriate recalibrating for D-aminoacid
580 read (ithep,*,end=111,err=111) nthetyp
581 do i=-nthetyp+1,nthetyp-1
582 read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
583 do j=0,nbend_kcc_Tb(i)
584 read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
589 & "Parameters of the valence-only potentials"
590 do i=-nthetyp+1,nthetyp-1
591 write (iout,'(2a)') "Type ",toronelet(i)
592 do k=0,nbend_kcc_Tb(i)
593 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
602 C write(iout,*) 'KURWA2'
605 C Read the parameters of the probability distribution/energy expression
606 C of the side chains.
609 cc write (iout,*) "tu dochodze",i
610 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
614 dsc_inv(i)=1.0D0/dsc(i)
625 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
626 & ((blower(k,l,1),l=1,k),k=1,3)
627 censc(1,1,-i)=censc(1,1,i)
628 censc(2,1,-i)=censc(2,1,i)
629 censc(3,1,-i)=-censc(3,1,i)
631 read (irotam,*,end=112,err=112) bsc(j,i)
632 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
633 & ((blower(k,l,j),l=1,k),k=1,3)
634 censc(1,j,-i)=censc(1,j,i)
635 censc(2,j,-i)=censc(2,j,i)
636 censc(3,j,-i)=-censc(3,j,i)
637 C BSC is amplitude of Gaussian
644 akl=akl+blower(k,m,j)*blower(l,m,j)
648 if (((k.eq.3).and.(l.ne.3))
649 & .or.((l.eq.3).and.(k.ne.3))) then
650 gaussc(k,l,j,-i)=-akl
651 gaussc(l,k,j,-i)=-akl
663 write (iout,'(/a)') 'Parameters of side-chain local geometry'
667 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
668 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
669 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
670 c write (iout,'(a,f10.4,4(16x,f10.4))')
671 c & 'Center ',(bsc(j,i),j=1,nlobi)
672 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
673 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
674 & 'log h',(bsc(j,i),j=1,nlobi)
675 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
676 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
683 c blower(k,l,j)=gaussc(ind,j,i)
688 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
689 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
696 C Read scrot parameters for potentials determined from all-atom AM1 calculations
697 C added by Urszula Kozlowska 07/11/2007
700 read (irotam,*,end=112,err=112)
702 read (irotam,*,end=112,err=112)
705 read(irotam,*,end=112,err=112) sc_parmin(j,i)
712 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
713 C interaction energy of the Gly, Ala, and Pro prototypes.
715 read (ifourier,*,end=115,err=115) nloctyp
716 SPLIT_FOURIERTOR = nloctyp.lt.0
717 nloctyp = iabs(nloctyp)
719 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
720 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
721 itype2loc(ntyp1)=nloctyp
722 iloctyp(nloctyp)=ntyp1
724 itype2loc(-i)=-itype2loc(i)
733 iloctyp(-i)=-iloctyp(i)
735 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
736 c write (iout,*) "nloctyp",nloctyp,
737 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
740 c write (iout,*) "NEWCORR",i
741 read (ifourier,*,end=115,err=115)
744 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
747 c write (iout,*) "NEWCORR BNEW1"
748 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
751 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
754 c write (iout,*) "NEWCORR BNEW2"
755 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
757 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
758 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
760 c write (iout,*) "NEWCORR CCNEW"
761 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
763 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
764 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
766 c write (iout,*) "NEWCORR DDNEW"
767 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
771 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
775 c write (iout,*) "NEWCORR EENEW1"
776 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
778 read (ifourier,*,end=115,err=115) e0new(ii,i)
780 c write (iout,*) (e0new(ii,i),ii=1,3)
782 c write (iout,*) "NEWCORR EENEW"
785 ccnew(ii,1,i)=ccnew(ii,1,i)/2
786 ccnew(ii,2,i)=ccnew(ii,2,i)/2
787 ddnew(ii,1,i)=ddnew(ii,1,i)/2
788 ddnew(ii,2,i)=ddnew(ii,2,i)/2
793 bnew1(ii,1,-i)= bnew1(ii,1,i)
794 bnew1(ii,2,-i)=-bnew1(ii,2,i)
795 bnew2(ii,1,-i)= bnew2(ii,1,i)
796 bnew2(ii,2,-i)=-bnew2(ii,2,i)
799 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
800 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
801 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
802 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
803 ccnew(ii,1,-i)=ccnew(ii,1,i)
804 ccnew(ii,2,-i)=-ccnew(ii,2,i)
805 ddnew(ii,1,-i)=ddnew(ii,1,i)
806 ddnew(ii,2,-i)=-ddnew(ii,2,i)
808 e0new(1,-i)= e0new(1,i)
809 e0new(2,-i)=-e0new(2,i)
810 e0new(3,-i)=-e0new(3,i)
812 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
813 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
814 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
815 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
819 write (iout,'(a)') "Coefficients of the multibody terms"
820 do i=-nloctyp+1,nloctyp-1
821 write (iout,*) "Type: ",onelet(iloctyp(i))
822 write (iout,*) "Coefficients of the expansion of B1"
824 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
826 write (iout,*) "Coefficients of the expansion of B2"
828 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
830 write (iout,*) "Coefficients of the expansion of C"
831 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
832 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
833 write (iout,*) "Coefficients of the expansion of D"
834 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
835 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
836 write (iout,*) "Coefficients of the expansion of E"
837 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
840 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
845 IF (SPLIT_FOURIERTOR) THEN
847 c write (iout,*) "NEWCORR TOR",i
848 read (ifourier,*,end=115,err=115)
851 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
854 c write (iout,*) "NEWCORR BNEW1 TOR"
855 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
858 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
861 c write (iout,*) "NEWCORR BNEW2 TOR"
862 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
864 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
865 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
867 c write (iout,*) "NEWCORR CCNEW TOR"
868 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
870 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
871 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
873 c write (iout,*) "NEWCORR DDNEW TOR"
874 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
878 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
882 c write (iout,*) "NEWCORR EENEW1 TOR"
883 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
885 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
887 c write (iout,*) (e0newtor(ii,i),ii=1,3)
889 c write (iout,*) "NEWCORR EENEW TOR"
892 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
893 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
894 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
895 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
900 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
901 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
902 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
903 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
906 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
907 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
908 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
909 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
911 e0newtor(1,-i)= e0newtor(1,i)
912 e0newtor(2,-i)=-e0newtor(2,i)
913 e0newtor(3,-i)=-e0newtor(3,i)
915 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
916 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
917 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
918 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
923 & "Single-body coefficients of the torsional potentials"
924 do i=-nloctyp+1,nloctyp-1
925 write (iout,*) "Type: ",onelet(iloctyp(i))
926 write (iout,*) "Coefficients of the expansion of B1tor"
928 write (iout,'(3hB1(,i1,1h),3f10.5)')
929 & j,(bnew1tor(k,j,i),k=1,3)
931 write (iout,*) "Coefficients of the expansion of B2tor"
933 write (iout,'(3hB2(,i1,1h),3f10.5)')
934 & j,(bnew2tor(k,j,i),k=1,3)
936 write (iout,*) "Coefficients of the expansion of Ctor"
937 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
938 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
939 write (iout,*) "Coefficients of the expansion of Dtor"
940 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
941 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
942 write (iout,*) "Coefficients of the expansion of Etor"
943 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
946 write (iout,'(1hE,2i1,2f10.5)')
947 & j,k,(eenewtor(l,j,k,i),l=1,2)
953 do i=-nloctyp+1,nloctyp-1
956 bnew1tor(ii,j,i)=bnew1(ii,j,i)
961 bnew2tor(ii,j,i)=bnew2(ii,j,i)
965 ccnewtor(ii,1,i)=ccnew(ii,1,i)
966 ccnewtor(ii,2,i)=ccnew(ii,2,i)
967 ddnewtor(ii,1,i)=ddnew(ii,1,i)
968 ddnewtor(ii,2,i)=ddnew(ii,2,i)
971 ENDIF !SPLIT_FOURIER_TOR
974 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
976 read (ifourier,*,end=115,err=115)
977 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
979 write (iout,*) 'Type ',onelet(iloctyp(i))
980 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
994 c B1tilde(1,i) = b(3)
995 c B1tilde(2,i) =-b(5)
996 c B1tilde(1,-i) =-b(3)
997 c B1tilde(2,-i) =b(5)
1004 cc B1tilde(1,i) = b(3,i)
1005 cc B1tilde(2,i) =-b(5,i)
1006 C B1tilde(1,-i) =-b(3,i)
1007 C B1tilde(2,-i) =b(5,i)
1008 cc b1tilde(1,i)=0.0d0
1009 cc b1tilde(2,i)=0.0d0
1017 CCold(1,1,i)= b(7,i)
1018 CCold(2,2,i)=-b(7,i)
1019 CCold(2,1,i)= b(9,i)
1020 CCold(1,2,i)= b(9,i)
1021 CCold(1,1,-i)= b(7,i)
1022 CCold(2,2,-i)=-b(7,i)
1023 CCold(2,1,-i)=-b(9,i)
1024 CCold(1,2,-i)=-b(9,i)
1029 c Ctilde(1,1,i)= CCold(1,1,i)
1030 c Ctilde(1,2,i)= CCold(1,2,i)
1031 c Ctilde(2,1,i)=-CCold(2,1,i)
1032 c Ctilde(2,2,i)=-CCold(2,2,i)
1034 c Ctilde(1,1,i)=0.0d0
1035 c Ctilde(1,2,i)=0.0d0
1036 c Ctilde(2,1,i)=0.0d0
1037 c Ctilde(2,2,i)=0.0d0
1038 DDold(1,1,i)= b(6,i)
1039 DDold(2,2,i)=-b(6,i)
1040 DDold(2,1,i)= b(8,i)
1041 DDold(1,2,i)= b(8,i)
1042 DDold(1,1,-i)= b(6,i)
1043 DDold(2,2,-i)=-b(6,i)
1044 DDold(2,1,-i)=-b(8,i)
1045 DDold(1,2,-i)=-b(8,i)
1050 c Dtilde(1,1,i)= DD(1,1,i)
1051 c Dtilde(1,2,i)= DD(1,2,i)
1052 c Dtilde(2,1,i)=-DD(2,1,i)
1053 c Dtilde(2,2,i)=-DD(2,2,i)
1055 c Dtilde(1,1,i)=0.0d0
1056 c Dtilde(1,2,i)=0.0d0
1057 c Dtilde(2,1,i)=0.0d0
1058 c Dtilde(2,2,i)=0.0d0
1059 EEold(1,1,i)= b(10,i)+b(11,i)
1060 EEold(2,2,i)=-b(10,i)+b(11,i)
1061 EEold(2,1,i)= b(12,i)-b(13,i)
1062 EEold(1,2,i)= b(12,i)+b(13,i)
1063 EEold(1,1,-i)= b(10,i)+b(11,i)
1064 EEold(2,2,-i)=-b(10,i)+b(11,i)
1065 EEold(2,1,-i)=-b(12,i)+b(13,i)
1066 EEold(1,2,-i)=-b(12,i)-b(13,i)
1067 c write(iout,*) "TU DOCHODZE"
1073 c ee(2,1,i)=ee(1,2,i)
1078 &"Coefficients of the cumulants (independent of valence angles)"
1079 do i=-nloctyp+1,nloctyp-1
1080 write (iout,*) 'Type ',onelet(iloctyp(i))
1082 write(iout,'(2f10.5)') B(3,i),B(5,i)
1084 write(iout,'(2f10.5)') B(2,i),B(4,i)
1087 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1091 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1095 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1100 C write (iout,*) 'KURWAKURWA'
1103 C Read torsional parameters in old format
1105 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1106 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1107 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1110 read (itorp,'(a)',end=113,err=113)
1112 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1118 write (iout,'(/a/)') 'Torsional constants:'
1121 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1122 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1128 C Read torsional parameters
1130 IF (TOR_MODE.eq.0) THEN
1132 read (itorp,*,end=113,err=113) ntortyp
1133 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1135 itype2loc(i)=itortyp(i)
1137 write (iout,*) 'ntortyp',ntortyp
1139 itype2loc(-i)=-itype2loc(i)
1141 itortyp(ntyp1)=ntortyp
1144 itortyp(i)=-itortyp(-i)
1146 c write (iout,*) 'ntortyp',ntortyp
1148 do j=-ntortyp+1,ntortyp-1
1149 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1151 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1152 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1155 do k=1,nterm(i,j,iblock)
1156 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1158 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1159 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1160 v0ij=v0ij+si*v1(k,i,j,iblock)
1163 do k=1,nlor(i,j,iblock)
1164 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1165 & vlor2(k,i,j),vlor3(k,i,j)
1166 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1169 v0(-i,-j,iblock)=v0ij
1175 write (iout,'(/a/)') 'Torsional constants:'
1178 write (iout,*) 'ityp',i,' jtyp',j
1179 write (iout,*) 'Fourier constants'
1180 do k=1,nterm(i,j,iblock)
1181 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1184 write (iout,*) 'Lorenz constants'
1185 do k=1,nlor(i,j,iblock)
1186 write (iout,'(3(1pe15.5))')
1187 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1193 C 6/23/01 Read parameters for double torsionals
1197 do j=-ntortyp+1,ntortyp-1
1198 do k=-ntortyp+1,ntortyp-1
1199 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1200 c write (iout,*) "OK onelett",
1203 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1204 & .or. t3.ne.toronelet(k)) then
1205 write (iout,*) "Error in double torsional parameter file",
1208 call MPI_Finalize(Ierror)
1210 stop "Error in double torsional parameter file"
1212 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1213 & ntermd_2(i,j,k,iblock)
1214 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1215 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1216 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1217 & ntermd_1(i,j,k,iblock))
1218 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1219 & ntermd_1(i,j,k,iblock))
1220 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1221 & ntermd_1(i,j,k,iblock))
1222 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1223 & ntermd_1(i,j,k,iblock))
1224 C Martix of D parameters for one dimesional foureir series
1225 do l=1,ntermd_1(i,j,k,iblock)
1226 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1227 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1228 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1229 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1230 c write(iout,*) "whcodze" ,
1231 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1233 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1234 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1235 & v2s(m,l,i,j,k,iblock),
1236 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1237 C Martix of D parameters for two dimesional fourier series
1238 do l=1,ntermd_2(i,j,k,iblock)
1240 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1241 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1242 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1243 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1252 write (iout,*) 'Constants for double torsionals'
1255 do j=-ntortyp+1,ntortyp-1
1256 do k=-ntortyp+1,ntortyp-1
1257 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1258 & ' nsingle',ntermd_1(i,j,k,iblock),
1259 & ' ndouble',ntermd_2(i,j,k,iblock)
1261 write (iout,*) 'Single angles:'
1262 do l=1,ntermd_1(i,j,k,iblock)
1263 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1264 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1265 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1266 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1269 write (iout,*) 'Pairs of angles:'
1270 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1271 do l=1,ntermd_2(i,j,k,iblock)
1272 write (iout,'(i5,20f10.5)')
1273 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1276 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1277 do l=1,ntermd_2(i,j,k,iblock)
1278 write (iout,'(i5,20f10.5)')
1279 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1280 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1289 ELSE IF (TOR_MODE.eq.1) THEN
1291 C read valence-torsional parameters
1292 read (itorp,*,end=113,err=113) ntortyp
1294 write (iout,*) "Valence-torsional parameters read in ntortyp",
1296 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1297 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1299 itortyp(i)=-itortyp(-i)
1301 do i=-ntortyp+1,ntortyp-1
1302 do j=-ntortyp+1,ntortyp-1
1303 C first we read the cos and sin gamma parameters
1304 read (itorp,'(13x,a)',end=113,err=113) string
1305 write (iout,*) i,j,string
1306 read (itorp,*,end=113,err=113)
1307 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1308 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1309 do k=1,nterm_kcc(j,i)
1310 do l=1,nterm_kcc_Tb(j,i)
1311 do ll=1,nterm_kcc_Tb(j,i)
1312 read (itorp,*,end=113,err=113) ii,jj,kk,
1313 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1321 c AL 4/8/16: Calculate coefficients from one-body parameters
1324 itortyp(i)=itype2loc(i)
1326 if (lprint) write (iout,*)
1327 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1329 do i=-ntortyp+1,ntortyp-1
1330 do j=-ntortyp+1,ntortyp-1
1333 do k=1,nterm_kcc_Tb(j,i)
1334 do l=1,nterm_kcc_Tb(j,i)
1335 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1336 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1337 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1338 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1341 do k=1,nterm_kcc_Tb(j,i)
1342 do l=1,nterm_kcc_Tb(j,i)
1344 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1345 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1346 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1347 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1349 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1350 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1351 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1352 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1356 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)
1362 if (tor_mode.gt.0 .and. lprint) then
1363 c Print valence-torsional parameters
1365 & "Parameters of the valence-torsional potentials"
1366 do i=-ntortyp+1,ntortyp-1
1367 do j=-ntortyp+1,ntortyp-1
1368 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1369 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1370 do k=1,nterm_kcc(j,i)
1371 do l=1,nterm_kcc_Tb(j,i)
1372 do ll=1,nterm_kcc_Tb(j,i)
1373 write (iout,'(3i5,2f15.4)')
1374 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1383 C Read of Side-chain backbone correlation parameters
1384 C Modified 11 May 2012 by Adasko
1387 read (isccor,*,end=119,err=119) nsccortyp
1388 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1390 isccortyp(i)=-isccortyp(-i)
1392 iscprol=isccortyp(20)
1393 c write (iout,*) 'ntortyp',ntortyp
1395 cc maxinter is maximum interaction sites
1399 read (isccor,*,end=119,err=119)
1400 &nterm_sccor(i,j),nlor_sccor(i,j)
1401 c write (iout,*) nterm_sccor(i,j)
1407 nterm_sccor(-i,j)=nterm_sccor(i,j)
1408 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1409 nterm_sccor(i,-j)=nterm_sccor(i,j)
1410 c write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1411 c & nterm_sccor(-i,-j),nterm_sccor(i,-j)
1412 do k=1,nterm_sccor(i,j)
1413 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1415 if (j.eq.iscprol) then
1416 if (i.eq.isccortyp(10)) then
1417 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1418 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1420 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1421 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1422 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1423 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1424 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1425 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1426 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1427 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1430 if (i.eq.isccortyp(10)) then
1431 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1432 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1434 if (j.eq.isccortyp(10)) then
1435 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1436 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1438 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1439 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)
1447 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1448 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1449 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1450 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1453 do k=1,nlor_sccor(i,j)
1454 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1455 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1456 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1457 &(1+vlor3sccor(k,i,j)**2)
1459 v0sccor(l,i,j)=v0ijsccor
1460 v0sccor(l,-i,j)=v0ijsccor1
1461 v0sccor(l,i,-j)=v0ijsccor2
1462 v0sccor(l,-i,-j)=v0ijsccor3
1468 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1472 write (iout,*) 'ityp',i,' jtyp',j
1473 write (iout,*) 'Fourier constants'
1474 do k=1,nterm_sccor(i,j)
1475 write (iout,'(2(1pe15.5))')
1476 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1478 write (iout,*) 'Lorenz constants'
1479 do k=1,nlor_sccor(i,j)
1480 write (iout,'(3(1pe15.5))')
1481 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1488 C Read electrostatic-interaction parameters
1491 write (iout,'(/a)') 'Electrostatic interaction constants:'
1492 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1493 & 'IT','JT','APP','BPP','AEL6','AEL3'
1495 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1496 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1497 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1498 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1503 app (i,j)=epp(i,j)*rri*rri
1504 bpp (i,j)=-2.0D0*epp(i,j)*rri
1505 ael6(i,j)=elpp6(i,j)*4.2D0**6
1506 ael3(i,j)=elpp3(i,j)*4.2D0**3
1507 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1508 & ael6(i,j),ael3(i,j)
1512 C Read side-chain interaction parameters.
1514 read (isidep,*,end=117,err=117) ipot,expon
1515 if (ipot.lt.1 .or. ipot.gt.5) then
1516 write (iout,'(2a)') 'Error while reading SC interaction',
1517 & 'potential file - unknown potential type.'
1521 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1522 & ', exponents are ',expon,2*expon
1523 goto (10,20,30,30,40) ipot
1524 C----------------------- LJ potential ---------------------------------
1525 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1526 & (sigma0(i),i=1,ntyp)
1528 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1529 write (iout,'(a/)') 'The epsilon array:'
1530 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1531 write (iout,'(/a)') 'One-body parameters:'
1532 write (iout,'(a,4x,a)') 'residue','sigma'
1533 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1536 C----------------------- LJK potential --------------------------------
1537 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1538 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1540 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1541 write (iout,'(a/)') 'The epsilon array:'
1542 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1543 write (iout,'(/a)') 'One-body parameters:'
1544 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1545 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1549 C---------------------- GB or BP potential -----------------------------
1551 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1553 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1554 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1555 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1556 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1558 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1559 C write(iout,*) "WARNING!!",i,ntyp
1560 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1562 C epslip(i,j)=epslip(i,j)+0.05d0
1565 C For the GB potential convert sigma'**2 into chi'
1568 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1572 write (iout,'(/a/)') 'Parameters of the BP potential:'
1573 write (iout,'(a/)') 'The epsilon array:'
1574 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1575 write (iout,'(/a)') 'One-body parameters:'
1576 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1578 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1579 & chip(i),alp(i),i=1,ntyp)
1582 C--------------------- GBV potential -----------------------------------
1583 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1584 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1585 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1587 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1588 write (iout,'(a/)') 'The epsilon array:'
1589 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1590 write (iout,'(/a)') 'One-body parameters:'
1591 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1592 & 's||/s_|_^2',' chip ',' alph '
1593 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1594 & sigii(i),chip(i),alp(i),i=1,ntyp)
1598 C-----------------------------------------------------------------------
1599 C Calculate the "working" parameters of SC interactions.
1603 epslip(i,j)=epslip(j,i)
1608 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1609 sigma(j,i)=sigma(i,j)
1610 rs0(i,j)=dwa16*sigma(i,j)
1614 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1615 & 'Working parameters of the SC interactions:',
1616 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1621 epsijlip=epslip(i,j)
1622 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1631 sigeps=dsign(1.0D0,epsij)
1633 aa_aq(i,j)=epsij*rrij*rrij
1634 bb_aq(i,j)=-sigeps*epsij*rrij
1635 aa_aq(j,i)=aa_aq(i,j)
1636 bb_aq(j,i)=bb_aq(i,j)
1637 sigeps=dsign(1.0D0,epsijlip)
1638 epsijlip=dabs(epsijlip)
1639 aa_lip(i,j)=epsijlip*rrij*rrij
1640 bb_lip(i,j)=-sigeps*epsijlip*rrij
1641 aa_lip(j,i)=aa_lip(i,j)
1642 bb_lip(j,i)=bb_lip(i,j)
1644 sigt1sq=sigma0(i)**2
1645 sigt2sq=sigma0(j)**2
1648 ratsig1=sigt2sq/sigt1sq
1649 ratsig2=1.0D0/ratsig1
1650 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1651 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1652 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1656 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1657 sigmaii(i,j)=rsum_max
1658 sigmaii(j,i)=rsum_max
1660 c sigmaii(i,j)=r0(i,j)
1661 c sigmaii(j,i)=r0(i,j)
1663 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1664 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1665 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1666 augm(i,j)=epsij*r_augm**(2*expon)
1667 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1674 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1675 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1676 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1681 C Define the SC-p interaction constants
1685 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1687 c aad(i,1)=0.3D0*4.0D0**12
1688 C Following line for constants currently implemented
1689 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1690 aad(i,1)=1.5D0*4.0D0**12
1691 c aad(i,1)=0.17D0*5.6D0**12
1693 C "Soft" SC-p repulsion
1695 C Following line for constants currently implemented
1696 c aad(i,1)=0.3D0*4.0D0**6
1697 C "Hard" SC-p repulsion
1698 bad(i,1)=3.0D0*4.0D0**6
1699 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1708 C 8/9/01 Read the SC-p interaction constants from file
1711 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1714 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1715 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1716 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1717 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1721 write (iout,'(/a)') "Parameters of SC-p interactions:"
1723 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1724 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1729 C Define the constants of the disulfide bridge
1733 c Old arbitrary potential - commented out.
1738 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1739 c energy surface of diethyl disulfide.
1740 c A. Liwo and U. Kozlowska, 11/24/03
1750 & write (iout,*) 'Dynamic formation/breaking of disulfides'
1752 ss_depth=ebr/wsc-0.25*eps(1,1)
1753 C write(iout,*) akcm,whpb,wsc,'KURWA'
1754 Ht=Ht/wsc-0.25*eps(1,1)
1763 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1768 write (iout,'(/a)') "Disulfide bridge parameters:"
1769 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1770 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1771 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1772 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1774 write (iout,'(a)') "Parameters of the 'trisulfide' potential"
1775 write(iout,*) "ATRISS",atriss
1776 write(iout,*) "BTRISS",btriss
1777 write(iout,*) "CTRISS",ctriss
1778 write(iout,*) "DTRISS",dtriss
1781 if (shield_mode.gt.0) then
1783 C VSolvSphere the volume of solving sphere
1785 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1786 C there will be no distinction between proline peptide group and normal peptide
1787 C group in case of shielding parameters
1788 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1789 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1790 write (iout,*) "VSolvSphere",VSolvSphere," VSolvSphere_div",
1792 C long axis of side chain
1794 long_r_sidechain(i)=vbldsc0(1,i)
1795 short_r_sidechain(i)=sigma0(i)
1800 111 write (iout,*) "Error reading bending energy parameters."
1802 112 write (iout,*) "Error reading rotamer energy parameters."
1804 113 write (iout,*) "Error reading torsional energy parameters."
1806 114 write (iout,*) "Error reading double torsional energy parameters."
1809 & "Error reading cumulant (multibody energy) parameters."
1811 116 write (iout,*) "Error reading electrostatic energy parameters."
1813 1161 write (iout,*) "Error reading lipid parameters."
1815 117 write (iout,*) "Error reading side chain interaction parameters."
1817 118 write (iout,*) "Error reading SCp interaction parameters."
1819 119 write (iout,*) "Error reading SCCOR parameters"
1821 121 write (iout,*) "Error reading bond parameters"
1824 call MPI_Finalize(Ierror)