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'
13 include 'COMMON.IOUNITS'
14 include 'COMMON.CHAIN'
15 include 'COMMON.INTERACT'
17 include 'COMMON.LOCAL'
18 include 'COMMON.TORSION'
19 include 'COMMON.FFIELD'
20 include 'COMMON.NAMES'
21 include 'COMMON.SBRIDGE'
22 include 'COMMON.WEIGHTS'
23 include 'COMMON.SCCOR'
24 include 'COMMON.SCROT'
25 include 'COMMON.SHIELD'
26 include 'COMMON.CONTROL'
27 include 'COMMON.DERIV'
28 include 'COMMON.OPTIM'
29 include 'COMMON.ENERGIES'
31 character*1 onelett(4) /"G","A","P","D"/
32 character*1 toronelet(-2:2) /"p","a","G","A","P"/
34 dimension blower(3,3,maxlob)
35 character*800 controlcard
36 character*256 bondname_t,thetname_t,rotname_t,torname_t,
37 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
43 double precision ip,mp
45 character*3 lancuch,ucase
49 if (me.eq.Master) then
50 call getenv("PRINT_PARM",lancuch)
51 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
55 C Set LPRINT=.TRUE. for debugging
56 dwa16=2.0d0**(1.0d0/6.0d0)
58 C Assign virtual-bond length
62 call card_concat(controlcard,.true.)
65 key = wname(i)(:ilen(wname(i)))
66 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
68 c print *,me," checkpoint 7"
70 if (me.eq.Master)write (iout,*) "iparm",iparm," myparm",myparm
72 c If reading not own parameters, skip assignment
73 call reada(controlcard,"D0CM",d0cm,3.78d0)
74 call reada(controlcard,"AKCM",akcm,15.1d0)
75 call reada(controlcard,"AKTH",akth,11.0d0)
76 call reada(controlcard,"AKCT",akct,12.0d0)
77 call reada(controlcard,"V1SS",v1ss,-1.08d0)
78 call reada(controlcard,"V2SS",v2ss,7.61d0)
79 call reada(controlcard,"V3SS",v3ss,13.7d0)
80 call reada(controlcard,"EBR",ebr,-5.50D0)
81 call reada(controlcard,"DTRISS",dtriss,1.0D0)
82 call reada(controlcard,"ATRISS",atriss,0.3D0)
83 call reada(controlcard,"BTRISS",btriss,0.02D0)
84 call reada(controlcard,"CTRISS",ctriss,1.0D0)
85 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
87 write(iout,*) "ATRISS",atriss
88 write(iout,*) "BTRISS",btriss
89 write(iout,*) "CTRISS",ctriss
90 write(iout,*) "DTRISS",dtriss
93 C dyn_ss_mask(i)=.false.
97 c Old arbitrary potential - commented out.
102 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
103 c energy surface of diethyl disulfide.
104 c A. Liwo and U. Kozlowska, 11/24/03
116 dyn_ssbond_ij(i,j)=1.0d300
119 call reada(controlcard,"HT",Ht,0.0D0)
121 C ss_depth=ebr/wsc-0.25*eps(1,1)
122 C write(iout,*) HT,wsc,eps(1,1),'KURWA'
123 C Ht=Ht/wsc-0.25*eps(1,1)
132 C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
135 if (iparm.eq.myparm) then
138 c Setup weights for UNRES
163 call card_concat(controlcard,.false.)
165 c Return if not own parameters
167 if (iparm.ne.myparm) return
169 call reads(controlcard,"BONDPAR",bondname_t,bondname)
170 open (ibond,file=bondname_t,status='old')
172 call reads(controlcard,"THETPAR",thetname_t,thetname)
173 open (ithep,file=thetname_t,status='old')
175 call reads(controlcard,"ROTPAR",rotname_t,rotname)
176 open (irotam,file=rotname_t,status='old')
178 call reads(controlcard,"TORPAR",torname_t,torname)
179 open (itorp,file=torname_t,status='old')
181 call reads(controlcard,"TORDPAR",tordname_t,tordname)
183 write (iout,*) "tor_mode",tor_mode
187 & open (itordp,file=tordname_t,status='old')
189 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
190 open (isccor,file=sccorname_t,status='old')
192 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
193 open (ifourier,file=fouriername_t,status='old')
195 call reads(controlcard,"ELEPAR",elename_t,elename)
196 open (ielep,file=elename_t,status='old')
198 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
199 open (isidep,file=sidename_t,status='old')
201 call reads(controlcard,"SCPPAR",scpname_t,scpname)
202 open (iscpp,file=scpname_t,status='old')
205 write (iout,*) "Parameter set:",iparm
206 write (iout,*) "Energy-term weights:"
208 write (iout,'(a16,f10.5)') wname(i),ww(i)
210 write (iout,*) "Sidechain potential file : ",
211 & sidename_t(:ilen(sidename_t))
213 write (iout,*) "SCp potential file : ",
214 & scpname_t(:ilen(scpname_t))
216 write (iout,*) "Electrostatic potential file : ",
217 & elename_t(:ilen(elename_t))
218 write (iout,*) "Cumulant coefficient file : ",
219 & fouriername_t(:ilen(fouriername_t))
220 write (iout,*) "Torsional parameter file : ",
221 & torname_t(:ilen(torname_t))
222 write (iout,*) "Double torsional parameter file : ",
223 & tordname_t(:ilen(tordname_t))
224 write (iout,*) "Backbone-rotamer parameter file : ",
225 & sccorname(:ilen(sccorname))
226 write (iout,*) "Bond & inertia constant file : ",
227 & bondname_t(:ilen(bondname_t))
228 write (iout,*) "Bending parameter file : ",
229 & thetname_t(:ilen(thetname_t))
230 write (iout,*) "Rotamer parameter file : ",
231 & rotname_t(:ilen(rotname_t))
234 c Read the virtual-bond parameters, masses, and moments of inertia
235 c and Stokes' radii of the peptide group and side chains
238 read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp
241 read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i)
242 dsc(i) = vbldsc0(1,i)
246 dsc_inv(i)=1.0D0/dsc(i)
250 read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk
252 read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
253 & aksc(j,i),abond0(j,i),j=1,nbondterm(i))
254 dsc(i) = vbldsc0(1,i)
258 dsc_inv(i)=1.0D0/dsc(i)
263 write(iout,'(/a/)')"Force constants virtual bonds:"
264 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
266 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
268 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
269 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
271 write (iout,'(13x,3f10.5)')
272 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
276 read(iliptranpar,*,end=1161,err=1161) pepliptran
278 read(iliptranpar,*,end=1161,err=1161) liptranene(i)
283 C Read the parameters of the probability distribution/energy expression
284 C of the virtual-bond valence angles theta
287 read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
288 & (bthet(j,i,1,1),j=1,2)
289 read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
290 read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
291 read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
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)
306 athet(1,i,-1,-1)=athet(1,-i,1,1)
307 athet(2,i,-1,-1)=-athet(2,-i,1,1)
308 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
309 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
310 athet(1,i,-1,1)=athet(1,-i,1,1)
311 athet(2,i,-1,1)=-athet(2,-i,1,1)
312 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
313 bthet(2,i,-1,1)=bthet(2,-i,1,1)
314 athet(1,i,1,-1)=-athet(1,-i,1,1)
315 athet(2,i,1,-1)=athet(2,-i,1,1)
316 bthet(1,i,1,-1)=bthet(1,-i,1,1)
317 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
322 polthet(j,i)=polthet(j,-i)
325 gthet(j,i)=gthet(j,-i)
331 c & 'Parameters of the virtual-bond valence angles:'
332 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
333 c & ' ATHETA0 ',' A1 ',' A2 ',
336 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
337 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
339 c write (iout,'(/a/9x,5a/79(1h-))')
340 c & 'Parameters of the expression for sigma(theta_c):',
341 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
342 c & ' ALPH3 ',' SIGMA0C '
344 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
345 c & (polthet(j,i),j=0,3),sigc0(i)
347 c write (iout,'(/a/9x,5a/79(1h-))')
348 c & 'Parameters of the second gaussian:',
349 c & ' THETA0 ',' SIGMA0 ',' G1 ',
352 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
353 c & sig0(i),(gthet(j,i),j=1,3)
356 & 'Parameters of the virtual-bond valence angles:'
357 write (iout,'(/a/9x,5a/79(1h-))')
358 & 'Coefficients of expansion',
359 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
360 & ' b1*10^1 ',' b2*10^1 '
362 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
363 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
364 & (10*bthet(j,i,1,1),j=1,2)
366 write (iout,'(/a/9x,5a/79(1h-))')
367 & 'Parameters of the expression for sigma(theta_c):',
368 & ' alpha0 ',' alph1 ',' alph2 ',
369 & ' alhp3 ',' sigma0c '
371 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
372 & (polthet(j,i),j=0,3),sigc0(i)
374 write (iout,'(/a/9x,5a/79(1h-))')
375 & 'Parameters of the second gaussian:',
376 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
379 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
380 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
384 IF (tor_mode.eq.0) THEN
386 C Read the parameters of Utheta determined from ab initio surfaces
387 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
389 read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
390 & ntheterm3,nsingle,ndouble
391 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
392 read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
394 ithetyp(i)=-ithetyp(-i)
397 do i=-maxthetyp,maxthetyp
398 do j=-maxthetyp,maxthetyp
399 do k=-maxthetyp,maxthetyp
400 aa0thet(i,j,k,iblock)=0.0d0
402 aathet(l,i,j,k,iblock)=0.0d0
406 bbthet(m,l,i,j,k,iblock)=0.0d0
407 ccthet(m,l,i,j,k,iblock)=0.0d0
408 ddthet(m,l,i,j,k,iblock)=0.0d0
409 eethet(m,l,i,j,k,iblock)=0.0d0
415 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
416 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
424 C write (iout,*) "KURWA1"
427 do j=-nthetyp,nthetyp
428 do k=-nthetyp,nthetyp
429 read (ithep,'(6a)',end=111,err=111) res1
430 c write(iout,*) res1,i,j,k
431 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
432 read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
434 read (ithep,*,end=111,err=111)
435 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
436 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
437 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
438 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
440 read (ithep,*,end=111,err=111)
441 & (((ffthet(llll,lll,ll,i,j,k,iblock),
442 & ffthet(lll,llll,ll,i,j,k,iblock),
443 & ggthet(llll,lll,ll,i,j,k,iblock)
444 & ,ggthet(lll,llll,ll,i,j,k,iblock),
445 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
449 C write(iout,*) "KURWA1.1"
451 C For dummy ends assign glycine-type coefficients of theta-only terms; the
452 C coefficients of theta-and-gamma-dependent terms are zero.
457 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
458 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
460 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
461 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
464 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
466 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
469 C write(iout,*) "KURWA1.5"
470 C Substitution for D aminoacids from symmetry.
473 do j=-nthetyp,nthetyp
474 do k=-nthetyp,nthetyp
475 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
477 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
481 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
482 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
483 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
484 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
490 ffthet(llll,lll,ll,i,j,k,iblock)=
491 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
492 ffthet(lll,llll,ll,i,j,k,iblock)=
493 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
494 ggthet(llll,lll,ll,i,j,k,iblock)=
495 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
496 ggthet(lll,llll,ll,i,j,k,iblock)=
497 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
507 C Control printout of the coefficients of virtual-bond-angle potentials
510 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
515 write (iout,'(//4a)')
516 & 'Type ',onelett(i),onelett(j),onelett(k)
517 write (iout,'(//a,10x,a)') " l","a[l]"
518 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
519 write (iout,'(i2,1pe15.5)')
520 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
522 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
523 & "b",l,"c",l,"d",l,"e",l
525 write (iout,'(i2,4(1pe15.5))') m,
526 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
527 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
531 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
532 & "f+",l,"f-",l,"g+",l,"g-",l
535 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
536 & ffthet(n,m,l,i,j,k,iblock),
537 & ffthet(m,n,l,i,j,k,iblock),
538 & ggthet(n,m,l,i,j,k,iblock),
539 & ggthet(m,n,l,i,j,k,iblock)
552 C here will be the apropriate recalibrating for D-aminoacid
553 read (ithep,*,end=111,err=111) nthetyp
554 do i=-nthetyp+1,nthetyp-1
555 read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
556 do j=0,nbend_kcc_Tb(i)
557 read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
562 & "Parameters of the valence-only potentials"
563 do i=-nthetyp+1,nthetyp-1
564 write (iout,'(2a)') "Type ",toronelet(i)
565 do k=0,nbend_kcc_Tb(i)
566 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
575 C write(iout,*) 'KURWA2'
578 C Read the parameters of the probability distribution/energy expression
579 C of the side chains.
582 cc write (iout,*) "tu dochodze",i
583 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
587 dsc_inv(i)=1.0D0/dsc(i)
598 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
599 & ((blower(k,l,1),l=1,k),k=1,3)
600 censc(1,1,-i)=censc(1,1,i)
601 censc(2,1,-i)=censc(2,1,i)
602 censc(3,1,-i)=-censc(3,1,i)
604 read (irotam,*,end=112,err=112) bsc(j,i)
605 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
606 & ((blower(k,l,j),l=1,k),k=1,3)
607 censc(1,j,-i)=censc(1,j,i)
608 censc(2,j,-i)=censc(2,j,i)
609 censc(3,j,-i)=-censc(3,j,i)
610 C BSC is amplitude of Gaussian
617 akl=akl+blower(k,m,j)*blower(l,m,j)
621 if (((k.eq.3).and.(l.ne.3))
622 & .or.((l.eq.3).and.(k.ne.3))) then
623 gaussc(k,l,j,-i)=-akl
624 gaussc(l,k,j,-i)=-akl
636 write (iout,'(/a)') 'Parameters of side-chain local geometry'
640 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
641 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
642 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
643 c write (iout,'(a,f10.4,4(16x,f10.4))')
644 c & 'Center ',(bsc(j,i),j=1,nlobi)
645 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
646 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
647 & 'log h',(bsc(j,i),j=1,nlobi)
648 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
649 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
656 c blower(k,l,j)=gaussc(ind,j,i)
661 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
662 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
669 C Read scrot parameters for potentials determined from all-atom AM1 calculations
670 C added by Urszula Kozlowska 07/11/2007
673 read (irotam,*,end=112,err=112)
675 read (irotam,*,end=112,err=112)
678 read(irotam,*,end=112,err=112) sc_parmin(j,i)
685 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
686 C interaction energy of the Gly, Ala, and Pro prototypes.
688 read (ifourier,*,end=115,err=115) nloctyp
690 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
691 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
692 itype2loc(ntyp1)=nloctyp
693 iloctyp(nloctyp)=ntyp1
695 itype2loc(-i)=-itype2loc(i)
699 itype2loc(i)=itortyp(i)
707 iloctyp(-i)=-iloctyp(i)
709 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
710 c write (iout,*) "nloctyp",nloctyp,
711 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
714 c write (iout,*) "NEWCORR",i
715 read (ifourier,*,end=115,err=115)
718 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
721 c write (iout,*) "NEWCORR BNEW1"
722 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
725 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
728 c write (iout,*) "NEWCORR BNEW2"
729 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
731 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
732 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
734 c write (iout,*) "NEWCORR CCNEW"
735 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
737 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
738 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
740 c write (iout,*) "NEWCORR DDNEW"
741 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
745 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
749 c write (iout,*) "NEWCORR EENEW1"
750 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
752 read (ifourier,*,end=115,err=115) e0new(ii,i)
754 c write (iout,*) (e0new(ii,i),ii=1,3)
756 c write (iout,*) "NEWCORR EENEW"
759 ccnew(ii,1,i)=ccnew(ii,1,i)/2
760 ccnew(ii,2,i)=ccnew(ii,2,i)/2
761 ddnew(ii,1,i)=ddnew(ii,1,i)/2
762 ddnew(ii,2,i)=ddnew(ii,2,i)/2
767 bnew1(ii,1,-i)= bnew1(ii,1,i)
768 bnew1(ii,2,-i)=-bnew1(ii,2,i)
769 bnew2(ii,1,-i)= bnew2(ii,1,i)
770 bnew2(ii,2,-i)=-bnew2(ii,2,i)
773 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
774 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
775 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
776 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
777 ccnew(ii,1,-i)=ccnew(ii,1,i)
778 ccnew(ii,2,-i)=-ccnew(ii,2,i)
779 ddnew(ii,1,-i)=ddnew(ii,1,i)
780 ddnew(ii,2,-i)=-ddnew(ii,2,i)
782 e0new(1,-i)= e0new(1,i)
783 e0new(2,-i)=-e0new(2,i)
784 e0new(3,-i)=-e0new(3,i)
786 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
787 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
788 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
789 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
793 write (iout,'(a)') "Coefficients of the multibody terms"
794 do i=-nloctyp+1,nloctyp-1
795 write (iout,*) "Type: ",onelet(iloctyp(i))
796 write (iout,*) "Coefficients of the expansion of B1"
798 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
800 write (iout,*) "Coefficients of the expansion of B2"
802 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
804 write (iout,*) "Coefficients of the expansion of C"
805 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
806 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
807 write (iout,*) "Coefficients of the expansion of D"
808 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
809 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
810 write (iout,*) "Coefficients of the expansion of E"
811 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
814 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
821 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
823 read (ifourier,*,end=115,err=115)
824 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
826 write (iout,*) 'Type ',onelet(iloctyp(i))
827 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
841 c B1tilde(1,i) = b(3)
842 c B1tilde(2,i) =-b(5)
843 c B1tilde(1,-i) =-b(3)
844 c B1tilde(2,-i) =b(5)
851 cc B1tilde(1,i) = b(3,i)
852 cc B1tilde(2,i) =-b(5,i)
853 C B1tilde(1,-i) =-b(3,i)
854 C B1tilde(2,-i) =b(5,i)
855 cc b1tilde(1,i)=0.0d0
856 cc b1tilde(2,i)=0.0d0
868 CCold(1,1,-i)= b(7,i)
869 CCold(2,2,-i)=-b(7,i)
870 CCold(2,1,-i)=-b(9,i)
871 CCold(1,2,-i)=-b(9,i)
876 c Ctilde(1,1,i)= CCold(1,1,i)
877 c Ctilde(1,2,i)= CCold(1,2,i)
878 c Ctilde(2,1,i)=-CCold(2,1,i)
879 c Ctilde(2,2,i)=-CCold(2,2,i)
881 c Ctilde(1,1,i)=0.0d0
882 c Ctilde(1,2,i)=0.0d0
883 c Ctilde(2,1,i)=0.0d0
884 c Ctilde(2,2,i)=0.0d0
889 DDold(1,1,-i)= b(6,i)
890 DDold(2,2,-i)=-b(6,i)
891 DDold(2,1,-i)=-b(8,i)
892 DDold(1,2,-i)=-b(8,i)
897 c Dtilde(1,1,i)= DD(1,1,i)
898 c Dtilde(1,2,i)= DD(1,2,i)
899 c Dtilde(2,1,i)=-DD(2,1,i)
900 c Dtilde(2,2,i)=-DD(2,2,i)
902 c Dtilde(1,1,i)=0.0d0
903 c Dtilde(1,2,i)=0.0d0
904 c Dtilde(2,1,i)=0.0d0
905 c Dtilde(2,2,i)=0.0d0
906 EEold(1,1,i)= b(10,i)+b(11,i)
907 EEold(2,2,i)=-b(10,i)+b(11,i)
908 EEold(2,1,i)= b(12,i)-b(13,i)
909 EEold(1,2,i)= b(12,i)+b(13,i)
910 EEold(1,1,-i)= b(10,i)+b(11,i)
911 EEold(2,2,-i)=-b(10,i)+b(11,i)
912 EEold(2,1,-i)=-b(12,i)+b(13,i)
913 EEold(1,2,-i)=-b(12,i)-b(13,i)
914 write(iout,*) "TU DOCHODZE"
920 c ee(2,1,i)=ee(1,2,i)
925 &"Coefficients of the cumulants (independent of valence angles)"
926 do i=-nloctyp+1,nloctyp-1
927 write (iout,*) 'Type ',onelet(iloctyp(i))
929 write(iout,'(2f10.5)') B(3,i),B(5,i)
931 write(iout,'(2f10.5)') B(2,i),B(4,i)
934 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
938 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
942 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
947 C write (iout,*) 'KURWAKURWA'
950 C Read torsional parameters in old format
952 read (itorp,*,end=113,err=113) ntortyp,nterm_old
953 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
954 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
957 read (itorp,'(a)',end=113,err=113)
959 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
965 write (iout,'(/a/)') 'Torsional constants:'
968 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
969 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
975 C Read torsional parameters
977 IF (TOR_MODE.eq.0) THEN
979 read (itorp,*,end=113,err=113) ntortyp
980 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
981 if (lprint)write (iout,*) 'ntortyp',ntortyp
984 itortyp(i)=-itortyp(-i)
986 c write (iout,*) 'ntortyp',ntortyp
988 do j=-ntortyp+1,ntortyp-1
989 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
991 nterm(-i,-j,iblock)=nterm(i,j,iblock)
992 nlor(-i,-j,iblock)=nlor(i,j,iblock)
995 do k=1,nterm(i,j,iblock)
996 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
998 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
999 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1000 v0ij=v0ij+si*v1(k,i,j,iblock)
1003 do k=1,nlor(i,j,iblock)
1004 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1005 & vlor2(k,i,j),vlor3(k,i,j)
1006 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1009 v0(-i,-j,iblock)=v0ij
1015 write (iout,'(/a/)') 'Torsional constants:'
1018 write (iout,*) 'ityp',i,' jtyp',j
1019 write (iout,*) 'Fourier constants'
1020 do k=1,nterm(i,j,iblock)
1021 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1024 write (iout,*) 'Lorenz constants'
1025 do k=1,nlor(i,j,iblock)
1026 write (iout,'(3(1pe15.5))')
1027 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1033 C 6/23/01 Read parameters for double torsionals
1037 do j=-ntortyp+1,ntortyp-1
1038 do k=-ntortyp+1,ntortyp-1
1039 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1040 c write (iout,*) "OK onelett",
1043 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1044 & .or. t3.ne.toronelet(k)) then
1046 & write (iout,*) "Error in double torsional parameter file",
1049 call MPI_Finalize(Ierror)
1051 stop "Error in double torsional parameter file"
1053 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1054 & ntermd_2(i,j,k,iblock)
1055 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1056 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1057 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1058 & ntermd_1(i,j,k,iblock))
1059 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1060 & ntermd_1(i,j,k,iblock))
1061 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1062 & ntermd_1(i,j,k,iblock))
1063 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1064 & ntermd_1(i,j,k,iblock))
1065 C Martix of D parameters for one dimesional foureir series
1066 do l=1,ntermd_1(i,j,k,iblock)
1067 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1068 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1069 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1070 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1071 c write(iout,*) "whcodze" ,
1072 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1074 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1075 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1076 & v2s(m,l,i,j,k,iblock),
1077 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1078 C Martix of D parameters for two dimesional fourier series
1079 do l=1,ntermd_2(i,j,k,iblock)
1081 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1082 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1083 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1084 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1093 write (iout,*) 'Constants for double torsionals'
1096 do j=-ntortyp+1,ntortyp-1
1097 do k=-ntortyp+1,ntortyp-1
1098 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1099 & ' nsingle',ntermd_1(i,j,k,iblock),
1100 & ' ndouble',ntermd_2(i,j,k,iblock)
1102 write (iout,*) 'Single angles:'
1103 do l=1,ntermd_1(i,j,k,iblock)
1104 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1105 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1106 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1107 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1110 write (iout,*) 'Pairs of angles:'
1111 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1112 do l=1,ntermd_2(i,j,k,iblock)
1113 write (iout,'(i5,20f10.5)')
1114 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1117 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1118 do l=1,ntermd_2(i,j,k,iblock)
1119 write (iout,'(i5,20f10.5)')
1120 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1121 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1130 ELSE IF (TOR_MODE.eq.1) THEN
1132 C read valence-torsional parameters
1133 read (itorp,*,end=113,err=113) ntortyp
1136 & write (iout,*) "Valence-torsional parameters read in ntortyp",
1138 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1140 & write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1142 itortyp(i)=-itortyp(-i)
1144 do i=-ntortyp+1,ntortyp-1
1145 do j=-ntortyp+1,ntortyp-1
1146 C first we read the cos and sin gamma parameters
1147 read (itorp,'(13x,a)',end=113,err=113) string
1148 if (me.eq.Master) write (iout,*) i,j,string
1149 read (itorp,*,end=113,err=113)
1150 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1151 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1152 do k=1,nterm_kcc(j,i)
1153 do l=1,nterm_kcc_Tb(j,i)
1154 do ll=1,nterm_kcc_Tb(j,i)
1155 read (itorp,*,end=113,err=113) ii,jj,kk,
1156 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1164 c AL 4/8/16: Calculate coefficients from one-body parameters
1167 itortyp(i)=itype2loc(i)
1170 if (me.eq.Master) write (iout,*)
1171 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1173 do i=-ntortyp+1,ntortyp-1
1174 do j=-ntortyp+1,ntortyp-1
1177 do k=1,nterm_kcc_Tb(j,i)
1178 do l=1,nterm_kcc_Tb(j,i)
1179 v1_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,1,j)
1180 & +bnew1(k,2,i)*bnew2(l,2,j)
1181 v2_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,2,j)
1182 & +bnew1(k,2,i)*bnew2(l,1,j)
1185 do k=1,nterm_kcc_Tb(j,i)
1186 do l=1,nterm_kcc_Tb(j,i)
1187 v1_kcc(k,l,2,i,j)=-0.25*(ccnew(k,1,i)*ddnew(l,1,j)
1188 & -ccnew(k,2,i)*ddnew(l,2,j))
1189 v2_kcc(k,l,2,i,j)=-0.25*(ccnew(k,2,i)*ddnew(l,1,j)
1190 & +ccnew(k,1,i)*ddnew(l,2,j))
1193 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)
1199 if (tor_mode.gt.0 .and. lprint) then
1200 c Print valence-torsional parameters
1202 & "Parameters of the valence-torsional potentials"
1203 do i=-ntortyp+1,ntortyp-1
1204 do j=-ntortyp+1,ntortyp-1
1205 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1206 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1207 do k=1,nterm_kcc(j,i)
1208 do l=1,nterm_kcc_Tb(j,i)
1209 do ll=1,nterm_kcc_Tb(j,i)
1210 write (iout,'(3i5,2f15.4)')
1211 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1220 C Read of Side-chain backbone correlation parameters
1221 C Modified 11 May 2012 by Adasko
1224 read (isccor,*,end=119,err=119) nsccortyp
1225 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1227 isccortyp(i)=-isccortyp(-i)
1229 iscprol=isccortyp(20)
1230 c write (iout,*) 'ntortyp',ntortyp
1232 cc maxinter is maximum interaction sites
1236 read (isccor,*,end=119,err=119)
1237 &nterm_sccor(i,j),nlor_sccor(i,j)
1238 c write (iout,*) nterm_sccor(i,j)
1244 nterm_sccor(-i,j)=nterm_sccor(i,j)
1245 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1246 nterm_sccor(i,-j)=nterm_sccor(i,j)
1247 c write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1248 c & nterm_sccor(-i,-j),nterm_sccor(i,-j)
1249 do k=1,nterm_sccor(i,j)
1250 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1252 if (j.eq.iscprol) then
1253 if (i.eq.isccortyp(10)) then
1254 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1255 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1257 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1258 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1259 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1260 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1261 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1262 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1263 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1264 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1267 if (i.eq.isccortyp(10)) then
1268 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1269 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1271 if (j.eq.isccortyp(10)) then
1272 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1273 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1275 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1276 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1277 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1278 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1279 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1280 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1284 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1285 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1286 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1287 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1290 do k=1,nlor_sccor(i,j)
1291 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1292 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1293 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1294 &(1+vlor3sccor(k,i,j)**2)
1296 v0sccor(l,i,j)=v0ijsccor
1297 v0sccor(l,-i,j)=v0ijsccor1
1298 v0sccor(l,i,-j)=v0ijsccor2
1299 v0sccor(l,-i,-j)=v0ijsccor3
1305 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1309 write (iout,*) 'ityp',i,' jtyp',j
1310 write (iout,*) 'Fourier constants'
1311 do k=1,nterm_sccor(i,j)
1312 write (iout,'(2(1pe15.5))')
1313 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1315 write (iout,*) 'Lorenz constants'
1316 do k=1,nlor_sccor(i,j)
1317 write (iout,'(3(1pe15.5))')
1318 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1325 C Read electrostatic-interaction parameters
1328 write (iout,'(/a)') 'Electrostatic interaction constants:'
1329 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1330 & 'IT','JT','APP','BPP','AEL6','AEL3'
1332 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1333 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1334 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1335 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1340 app (i,j)=epp(i,j)*rri*rri
1341 bpp (i,j)=-2.0D0*epp(i,j)*rri
1342 ael6(i,j)=elpp6(i,j)*4.2D0**6
1343 ael3(i,j)=elpp3(i,j)*4.2D0**3
1344 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1345 & ael6(i,j),ael3(i,j)
1349 C Read side-chain interaction parameters.
1351 read (isidep,*) ipot,expon
1352 if (ipot.lt.1 .or. ipot.gt.6) then
1354 & write (iout,'(2a)') 'Error while reading SC interaction',
1355 & 'potential file - unknown potential type.'
1360 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1361 & ', exponents are ',expon,2*expon
1362 goto (10,20,30,30,40,50,60) ipot
1363 C----------------------- LJ potential ---------------------------------
1364 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1366 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1367 write (iout,'(a/)') 'The epsilon array:'
1368 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1369 write (iout,'(/a)') 'One-body parameters:'
1370 write (iout,'(a,4x,a)') 'residue','sigma'
1371 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1374 C----------------------- LJK potential --------------------------------
1375 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1376 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1378 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1379 write (iout,'(a/)') 'The epsilon array:'
1380 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1381 write (iout,'(/a)') 'One-body parameters:'
1382 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1383 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1387 C---------------------- GB or BP potential -----------------------------
1388 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1389 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
1391 C For the GB potential convert sigma'**2 into chi'
1394 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1398 write (iout,'(/a/)') 'Parameters of the BP potential:'
1399 write (iout,'(a/)') 'The epsilon array:'
1400 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1401 write (iout,'(/a)') 'One-body parameters:'
1402 write (iout,'(a,4x,5a)') 'residue',' sigma ','s||/s_|_^2',
1403 & ' chip0 ',' chip ',' alph '
1404 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),sigii(i),
1405 & chip0(i),chip(i),alp(i),i=1,ntyp)
1408 C--------------------- GBV potential -----------------------------------
1409 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1410 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1411 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1413 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1414 write (iout,'(a/)') 'The epsilon array:'
1415 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1416 write (iout,'(/a)') 'One-body parameters:'
1417 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1418 & 's||/s_|_^2',' chip ',' alph '
1419 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1420 & sigii(i),chip(i),alp(i),i=1,ntyp)
1423 C--------------------- Momo potential -----------------------------------
1425 read (isidep,*) (icharge(i),i=1,ntyp)
1426 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
1429 c! write (*,*) "Im in ", i, " ", j
1431 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
1432 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
1433 & chis(i,j),chis(j,i),
1434 & nstate(i,j),(wstate(k,i,j),k=1,4),
1439 & dtail(1,i,j),dtail(2,i,j),
1440 & epshead(i,j),sig0head(i,j),
1441 & rborn(i,j),rborn(j,i),
1442 & (wqdip(k,i,j),k=1,2),wquad(i,j),
1443 & alphapol(i,j),alphapol(j,i),
1444 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
1447 c! write (*,*) "nstate gly-gly", nstate(10,10)
1448 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
1449 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
1450 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
1453 c! DISABLE IT AT >>YOUR OWN PERIL<<
1458 sigma(i,j) = sigma(j,i)
1459 nstate(i,j) = nstate(j,i)
1460 sigmap1(i,j) = sigmap1(j,i)
1461 sigmap2(i,j) = sigmap2(j,i)
1462 sigiso1(i,j) = sigiso1(j,i)
1463 sigiso2(i,j) = sigiso2(j,i)
1466 alphasur(k,i,j) = alphasur(k,j,i)
1467 wstate(k,i,j) = wstate(k,j,i)
1468 alphiso(k,i,j) = alphiso(k,j,i)
1471 dhead(2,1,i,j) = dhead(1,1,j,i)
1472 dhead(2,2,i,j) = dhead(1,2,j,i)
1473 dhead(1,1,i,j) = dhead(2,1,j,i)
1474 dhead(1,2,i,j) = dhead(2,2,j,i)
1475 dtail(1,i,j) = dtail(1,j,i)
1476 dtail(2,i,j) = dtail(2,j,i)
1479 c! write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
1480 c! write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
1481 c! dhead(l,k,i,j) = dhead(k,l,j,i)
1485 epshead(i,j) = epshead(j,i)
1486 sig0head(i,j) = sig0head(j,i)
1489 wqdip(k,i,j) = wqdip(k,j,i)
1492 wquad(i,j) = wquad(j,i)
1493 epsintab(i,j) = epsintab(j,i)
1498 if (.not.lprint) goto 70
1500 & "Parameters of the new physics-based SC-SC interaction potential"
1501 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
1502 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
1505 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
1506 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
1507 & chipp(i,j),chipp(j,i)
1510 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
1511 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
1515 write (iout,'(2(a3,1x),8(0pf10.3))')
1516 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
1517 & sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
1520 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
1521 & ' wst2',' wst3',' wst4',' dh11',' dh21',
1522 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
1526 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
1527 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
1528 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1529 & epshead(i,j),sig0head(i,j)
1532 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
1533 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
1537 write (iout,'(2(a3,1x),2i4,5f10.3)')
1538 & restyp(i),restyp(j),icharge(i),icharge(j),
1539 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
1542 write (iout,'(/12a)') 'Residues',
1544 & ' alphpol2',' alphiso1',' alpiso2',
1545 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
1549 write (iout,'(2(a3,1x),11f10.3)')
1550 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
1551 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
1556 C For the GB potential convert sigma'**2 into chi'
1558 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1561 write (iout,'(/a/)') 'Parameters of the BP potential:'
1562 write (iout,'(a/)') 'The epsilon array:'
1563 CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1564 write (iout,'(/a)') 'One-body parameters:'
1565 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1567 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1568 & chip(i),alp(i),i=1,ntyp)
1572 C-----------------------------------------------------------------------
1573 C Calculate the "working" parameters of SC interactions.
1578 alphasur(k,j,i)=alphasur(k,i,j)
1579 alphiso(k,j,i)=alphiso(k,i,j)
1580 wstate(k,j,i)=wstate(k,i,j)
1583 wqdip(k,j,i)=wqdip(k,i,j)
1587 dhead(l,k,j,i)=dhead(l,k,i,j)
1595 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1596 sigma(j,i)=sigma(i,j)
1597 rs0(i,j)=dwa16*sigma(i,j)
1604 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1605 & 'Working parameters of the SC interactions:',
1606 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1611 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 )
1621 sigeps=dsign(1.0D0,epsij)
1623 aa(i,j)=epsij*rrij*rrij
1624 bb(i,j)=-sigeps*epsij*rrij
1627 IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
1628 sigt1sq=sigma0(i)**2
1629 sigt2sq=sigma0(j)**2
1632 ratsig1=sigt2sq/sigt1sq
1633 ratsig2=1.0D0/ratsig1
1634 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1635 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1636 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1640 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1641 sigmaii(i,j)=rsum_max
1642 sigmaii(j,i)=rsum_max
1644 c sigmaii(i,j)=r0(i,j)
1645 c sigmaii(j,i)=r0(i,j)
1647 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1648 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1649 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1650 augm(i,j)=epsij*r_augm**(2*expon)
1651 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1659 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1660 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1661 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1663 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1665 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1666 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1667 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1668 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1669 & chis(i,j),chis(j,i),
1670 & nstate(i,j),(wstate(k,i,j),k=1,4),
1671 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1672 & epshead(i,j),sig0head(i,j),
1673 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1674 & alphapol(i,j),alphapol(j,i),
1675 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1682 C Define the SC-p interaction constants
1686 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1688 c aad(i,1)=0.3D0*4.0D0**12
1689 C Following line for constants currently implemented
1690 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1691 aad(i,1)=1.5D0*4.0D0**12
1692 c aad(i,1)=0.17D0*5.6D0**12
1694 C "Soft" SC-p repulsion
1696 C Following line for constants currently implemented
1697 c aad(i,1)=0.3D0*4.0D0**6
1698 C "Hard" SC-p repulsion
1699 bad(i,1)=3.0D0*4.0D0**6
1700 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1709 C 8/9/01 Read the SC-p interaction constants from file
1712 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1715 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1716 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1717 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1718 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1722 write (iout,*) "Parameters of SC-p interactions:"
1724 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1725 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1730 C Define the constants of the disulfide bridge
1734 c Old arbitrary potential - commented out.
1739 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1740 c energy surface of diethyl disulfide.
1741 c A. Liwo and U. Kozlowska, 11/24/03
1750 if (lprint) write (iout,*) dyn_ss,'dyndyn'
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
1767 if (me.eq.Master) then
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,
1776 if (shield_mode.gt.0) then
1778 C VSolvSphere the volume of solving sphere
1780 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1781 C there will be no distinction between proline peptide group and normal peptide
1782 C group in case of shielding parameters
1783 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1784 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1785 if (me.eq.Master) write (iout,*) VSolvSphere,VSolvSphere_div
1786 C long axis of side chain
1788 long_r_sidechain(i)=vbldsc0(1,i)
1789 short_r_sidechain(i)=sigma0(i)
1794 111 if (mq.eq.Masster)
1795 & write (iout,*) "Error reading bending energy parameters."
1797 112 if (me.eq.Master)
1798 & write (iout,*) "Error reading rotamer energy parameters."
1800 113 if (me.eq.Master)
1801 & write (iout,*) "Error reading torsional energy parameters."
1803 114 if (me.eq.Master)
1804 &write (iout,*) "Error reading double torsional energy parameters."
1806 115 if (me.eq.Master) write (iout,*)
1807 & "Error reading cumulant (multibody energy) parameters."
1809 116 if (mw.eq.Master)
1810 & write (iout,*) "Error reading electrostatic energy parameters."
1812 1161 if (me.eq.Master)
1813 & write (iout,*) "Error reading lipid energy parameters."
1815 117 if (me.eq.Master)
1816 & write (iout,*) "Error reading side chain interaction parameters."
1818 118 if (me.eq.Master)
1819 & write (iout,*) "Error reading SCp interaction parameters."
1821 119 if (me.eq.Master)
1822 & write (iout,*) "Error reading SCCOR parameters"
1824 121 if (me.eq.Master)
1825 & write (iout,*) "Error reading bond parameters"
1828 call MPI_Finalize(Ierror)