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'
27 character*1 onelett(4) /"G","A","P","D"/
28 character*1 toronelet(-2:2) /"p","a","G","A","P"/
30 dimension blower(3,3,maxlob)
31 character*900 controlcard
32 character*256 bondname_t,thetname_t,rotname_t,torname_t,
33 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
34 & sccorname_t,tubename_t
39 double precision ip,mp
41 C write (iout,*) "KURWA"
45 C Set LPRINT=.TRUE. for debugging
46 dwa16=2.0d0**(1.0d0/6.0d0)
49 C Assign virtual-bond length
53 call card_concat(controlcard,.true.)
56 key = wname(i)(:ilen(wname(i)))
57 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
60 write (iout,*) "iparm",iparm," myparm",myparm
61 c If reading not own parameters, skip assignment
62 call reada(controlcard,"D0CM",d0cm,3.78d0)
63 call reada(controlcard,"AKCM",akcm,15.1d0)
64 call reada(controlcard,"AKTH",akth,11.0d0)
65 call reada(controlcard,"AKCT",akct,12.0d0)
66 call reada(controlcard,"V1SS",v1ss,-1.08d0)
67 call reada(controlcard,"V2SS",v2ss,7.61d0)
68 call reada(controlcard,"V3SS",v3ss,13.7d0)
69 call reada(controlcard,"EBR",ebr,-5.50D0)
70 call reada(controlcard,"DTRISS",dtriss,1.0D0)
71 call reada(controlcard,"ATRISS",atriss,0.3D0)
72 call reada(controlcard,"BTRISS",btriss,0.02D0)
73 call reada(controlcard,"CTRISS",ctriss,1.0D0)
74 call reada(controlcard,"LIPSCALE",lipscale,1.3D0)
75 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
76 write(iout,*) "ATRISS",atriss
77 write(iout,*) "BTRISS",btriss
78 write(iout,*) "CTRISS",ctriss
79 write(iout,*) "DTRISS",dtriss
82 C dyn_ss_mask(i)=.false.
86 c Old arbitrary potential - commented out.
91 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
92 c energy surface of diethyl disulfide.
93 c A. Liwo and U. Kozlowska, 11/24/03
105 dyn_ssbond_ij(i,j)=1.0d300
108 call reada(controlcard,"HT",Ht,0.0D0)
110 C ss_depth=ebr/wsc-0.25*eps(1,1)
111 C write(iout,*) HT,wsc,eps(1,1),'KURWA'
112 C Ht=Ht/wsc-0.25*eps(1,1)
121 C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
124 if (iparm.eq.myparm .or. .not.separate_parset) then
127 c Setup weights for UNRES
153 call card_concat(controlcard,.false.)
155 c Return if not own parameters
157 if (iparm.ne.myparm .and. separate_parset) return
159 call reads(controlcard,"BONDPAR",bondname_t,bondname)
160 open (ibond,file=bondname_t,status='old')
162 call reads(controlcard,"THETPAR",thetname_t,thetname)
163 open (ithep,file=thetname_t,status='old')
165 call reads(controlcard,"ROTPAR",rotname_t,rotname)
166 open (irotam,file=rotname_t,status='old')
168 call reads(controlcard,"TORPAR",torname_t,torname)
169 open (itorp,file=torname_t,status='old')
171 call reads(controlcard,"TORDPAR",tordname_t,tordname)
172 open (itordp,file=tordname_t,status='old')
174 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
175 open (isccor,file=sccorname_t,status='old')
177 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
178 open (ifourier,file=fouriername_t,status='old')
180 call reads(controlcard,"ELEPAR",elename_t,elename)
181 open (ielep,file=elename_t,status='old')
183 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
184 open (isidep,file=sidename_t,status='old')
186 call reads(controlcard,"SCPPAR",scpname_t,scpname)
187 open (iscpp,file=scpname_t,status='old')
189 call reads(controlcard,"TUBEPAR",tubename_t,tubename)
190 write(iout,*) tubename_t
191 write(iout,*) tubename
192 open (itube,file=tubename_t,status='old')
195 write (iout,*) "Parameter set:",iparm
196 write (iout,*) "Energy-term weights:"
198 write (iout,'(a16,f10.5)') wname(i),ww(i)
200 write (iout,*) "Sidechain potential file : ",
201 & sidename_t(:ilen(sidename_t))
203 write (iout,*) "SCp potential file : ",
204 & scpname_t(:ilen(scpname_t))
206 write (iout,*) "Electrostatic potential file : ",
207 & elename_t(:ilen(elename_t))
208 write (iout,*) "Cumulant coefficient file : ",
209 & fouriername_t(:ilen(fouriername_t))
210 write (iout,*) "Torsional parameter file : ",
211 & torname_t(:ilen(torname_t))
212 write (iout,*) "Double torsional parameter file : ",
213 & tordname_t(:ilen(tordname_t))
214 write (iout,*) "Backbone-rotamer parameter file : ",
215 & sccorname(:ilen(sccorname))
216 write (iout,*) "Bond & inertia constant file : ",
217 & bondname_t(:ilen(bondname_t))
218 write (iout,*) "Bending parameter file : ",
219 & thetname_t(:ilen(thetname_t))
220 write (iout,*) "Rotamer parameter file : ",
221 & rotname_t(:ilen(rotname_t))
224 c Read the virtual-bond parameters, masses, and moments of inertia
225 c and Stokes' radii of the peptide group and side chains
228 read (ibond,*) vbldp0,vbldpdum,akp
231 read (ibond,*) vbldsc0(1,i),aksc(1,i)
232 dsc(i) = vbldsc0(1,i)
236 dsc_inv(i)=1.0D0/dsc(i)
240 read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk
242 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
244 dsc(i) = vbldsc0(1,i)
248 dsc_inv(i)=1.0D0/dsc(i)
253 write(iout,'(/a/)')"Force constants virtual bonds:"
254 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
256 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
258 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
259 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
261 write (iout,'(13x,3f10.5)')
262 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
266 read(iliptranpar,*) pepliptran
268 read(iliptranpar,*) liptranene(i)
273 C Read the parameters of the probability distribution/energy expression
274 C of the virtual-bond valence angles theta
277 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
278 & (bthet(j,i,1,1),j=1,2)
279 read (ithep,*) (polthet(j,i),j=0,3)
280 read (ithep,*) (gthet(j,i),j=1,3)
281 read (ithep,*) theta0(i),sig0(i),sigc0(i)
285 athet(1,i,1,-1)=athet(1,i,1,1)
286 athet(2,i,1,-1)=athet(2,i,1,1)
287 bthet(1,i,1,-1)=-bthet(1,i,1,1)
288 bthet(2,i,1,-1)=-bthet(2,i,1,1)
289 athet(1,i,-1,1)=-athet(1,i,1,1)
290 athet(2,i,-1,1)=-athet(2,i,1,1)
291 bthet(1,i,-1,1)=bthet(1,i,1,1)
292 bthet(2,i,-1,1)=bthet(2,i,1,1)
296 athet(1,i,-1,-1)=athet(1,-i,1,1)
297 athet(2,i,-1,-1)=-athet(2,-i,1,1)
298 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
299 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
300 athet(1,i,-1,1)=athet(1,-i,1,1)
301 athet(2,i,-1,1)=-athet(2,-i,1,1)
302 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
303 bthet(2,i,-1,1)=bthet(2,-i,1,1)
304 athet(1,i,1,-1)=-athet(1,-i,1,1)
305 athet(2,i,1,-1)=athet(2,-i,1,1)
306 bthet(1,i,1,-1)=bthet(1,-i,1,1)
307 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
312 polthet(j,i)=polthet(j,-i)
315 gthet(j,i)=gthet(j,-i)
321 c & 'Parameters of the virtual-bond valence angles:'
322 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
323 c & ' ATHETA0 ',' A1 ',' A2 ',
326 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
327 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
329 c write (iout,'(/a/9x,5a/79(1h-))')
330 c & 'Parameters of the expression for sigma(theta_c):',
331 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
332 c & ' ALPH3 ',' SIGMA0C '
334 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
335 c & (polthet(j,i),j=0,3),sigc0(i)
337 c write (iout,'(/a/9x,5a/79(1h-))')
338 c & 'Parameters of the second gaussian:',
339 c & ' THETA0 ',' SIGMA0 ',' G1 ',
342 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
343 c & sig0(i),(gthet(j,i),j=1,3)
346 & 'Parameters of the virtual-bond valence angles:'
347 write (iout,'(/a/9x,5a/79(1h-))')
348 & 'Coefficients of expansion',
349 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
350 & ' b1*10^1 ',' b2*10^1 '
352 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
353 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
354 & (10*bthet(j,i,1,1),j=1,2)
356 write (iout,'(/a/9x,5a/79(1h-))')
357 & 'Parameters of the expression for sigma(theta_c):',
358 & ' alpha0 ',' alph1 ',' alph2 ',
359 & ' alhp3 ',' sigma0c '
361 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
362 & (polthet(j,i),j=0,3),sigc0(i)
364 write (iout,'(/a/9x,5a/79(1h-))')
365 & 'Parameters of the second gaussian:',
366 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
369 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
370 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
375 C Read the parameters of Utheta determined from ab initio surfaces
376 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
378 write (iout,*) "tu dochodze"
379 read (ithep,*) nthetyp,ntheterm,ntheterm2,
380 & ntheterm3,nsingle,ndouble
381 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
382 read (ithep,*) (ithetyp(i),i=1,ntyp1)
384 ithetyp(i)=-ithetyp(-i)
386 write (iout,*) "tu dochodze"
388 do i=-maxthetyp,maxthetyp
389 do j=-maxthetyp,maxthetyp
390 do k=-maxthetyp,maxthetyp
391 aa0thet(i,j,k,iblock)=0.0d0
393 aathet(l,i,j,k,iblock)=0.0d0
397 bbthet(m,l,i,j,k,iblock)=0.0d0
398 ccthet(m,l,i,j,k,iblock)=0.0d0
399 ddthet(m,l,i,j,k,iblock)=0.0d0
400 eethet(m,l,i,j,k,iblock)=0.0d0
406 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
407 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
415 C write (iout,*) "KURWA1"
418 do j=-nthetyp,nthetyp
419 do k=-nthetyp,nthetyp
420 read (ithep,'(6a)') res1
421 write(iout,*) res1,i,j,k
422 read (ithep,*) aa0thet(i,j,k,iblock)
423 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
425 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
426 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
427 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
428 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
431 & (((ffthet(llll,lll,ll,i,j,k,iblock),
432 & ffthet(lll,llll,ll,i,j,k,iblock),
433 & ggthet(llll,lll,ll,i,j,k,iblock)
434 & ,ggthet(lll,llll,ll,i,j,k,iblock),
435 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
439 C write(iout,*) "KURWA1.1"
441 C For dummy ends assign glycine-type coefficients of theta-only terms; the
442 C coefficients of theta-and-gamma-dependent terms are zero.
447 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
448 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
450 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
451 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
454 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
456 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
459 C write(iout,*) "KURWA1.5"
460 C Substitution for D aminoacids from symmetry.
463 do j=-nthetyp,nthetyp
464 do k=-nthetyp,nthetyp
465 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
467 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
471 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
472 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
473 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
474 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
480 ffthet(llll,lll,ll,i,j,k,iblock)=
481 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
482 ffthet(lll,llll,ll,i,j,k,iblock)=
483 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
484 ggthet(llll,lll,ll,i,j,k,iblock)=
485 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
486 ggthet(lll,llll,ll,i,j,k,iblock)=
487 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
497 C Control printout of the coefficients of virtual-bond-angle potentials
500 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
504 write (iout,'(//4a)')
505 & 'Type ',onelett(i),onelett(j),onelett(k)
506 write (iout,'(//a,10x,a)') " l","a[l]"
507 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
508 write (iout,'(i2,1pe15.5)')
509 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
511 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
512 & "b",l,"c",l,"d",l,"e",l
514 write (iout,'(i2,4(1pe15.5))') m,
515 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
516 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
520 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
521 & "f+",l,"f-",l,"g+",l,"g-",l
524 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
525 & ffthet(n,m,l,i,j,k,iblock),
526 & ffthet(m,n,l,i,j,k,iblock),
527 & ggthet(n,m,l,i,j,k,iblock),
528 & ggthet(m,n,l,i,j,k,iblock)
538 C write(iout,*) 'KURWA2'
541 C Read the parameters of the probability distribution/energy expression
542 C of the side chains.
545 cc write (iout,*) "tu dochodze",i
546 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
550 dsc_inv(i)=1.0D0/dsc(i)
561 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
562 censc(1,1,-i)=censc(1,1,i)
563 censc(2,1,-i)=censc(2,1,i)
564 censc(3,1,-i)=-censc(3,1,i)
566 read (irotam,*) bsc(j,i)
567 read (irotam,*) (censc(k,j,i),k=1,3),
568 & ((blower(k,l,j),l=1,k),k=1,3)
569 censc(1,j,-i)=censc(1,j,i)
570 censc(2,j,-i)=censc(2,j,i)
571 censc(3,j,-i)=-censc(3,j,i)
572 C BSC is amplitude of Gaussian
579 akl=akl+blower(k,m,j)*blower(l,m,j)
583 if (((k.eq.3).and.(l.ne.3))
584 & .or.((l.eq.3).and.(k.ne.3))) then
585 gaussc(k,l,j,-i)=-akl
586 gaussc(l,k,j,-i)=-akl
598 write (iout,'(/a)') 'Parameters of side-chain local geometry'
602 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
603 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
604 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
605 c write (iout,'(a,f10.4,4(16x,f10.4))')
606 c & 'Center ',(bsc(j,i),j=1,nlobi)
607 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
608 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
609 & 'log h',(bsc(j,i),j=1,nlobi)
610 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
611 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
618 c blower(k,l,j)=gaussc(ind,j,i)
623 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
624 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
631 C Read scrot parameters for potentials determined from all-atom AM1 calculations
632 C added by Urszula Kozlowska 07/11/2007
640 read(irotam,*) sc_parmin(j,i)
646 C write (iout,*) 'KURWAKURWA'
649 C Read torsional parameters in old format
651 read (itorp,*) ntortyp,nterm_old
652 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
653 read (itorp,*) (itortyp(i),i=1,ntyp)
658 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
664 write (iout,'(/a/)') 'Torsional constants:'
667 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
668 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
676 C Read torsional parameters
678 read (itorp,*) ntortyp
679 read (itorp,*) (itortyp(i),i=1,ntyp)
680 write (iout,*) 'ntortyp',ntortyp
683 itortyp(i)=-itortyp(-i)
685 c write (iout,*) 'ntortyp',ntortyp
687 do j=-ntortyp+1,ntortyp-1
688 read (itorp,*) nterm(i,j,iblock),
690 nterm(-i,-j,iblock)=nterm(i,j,iblock)
691 nlor(-i,-j,iblock)=nlor(i,j,iblock)
694 do k=1,nterm(i,j,iblock)
695 read (itorp,*) kk,v1(k,i,j,iblock),
697 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
698 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
699 v0ij=v0ij+si*v1(k,i,j,iblock)
702 do k=1,nlor(i,j,iblock)
703 read (itorp,*) kk,vlor1(k,i,j),
704 & vlor2(k,i,j),vlor3(k,i,j)
705 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
708 v0(-i,-j,iblock)=v0ij
714 write (iout,'(/a/)') 'Torsional constants:'
717 write (iout,*) 'ityp',i,' jtyp',j
718 write (iout,*) 'Fourier constants'
719 do k=1,nterm(i,j,iblock)
720 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
723 write (iout,*) 'Lorenz constants'
724 do k=1,nlor(i,j,iblock)
725 write (iout,'(3(1pe15.5))')
726 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
732 C 6/23/01 Read parameters for double torsionals
736 do j=-ntortyp+1,ntortyp-1
737 do k=-ntortyp+1,ntortyp-1
738 read (itordp,'(3a1)') t1,t2,t3
739 c write (iout,*) "OK onelett",
742 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
743 & .or. t3.ne.toronelet(k)) then
744 write (iout,*) "Error in double torsional parameter file",
747 call MPI_Finalize(Ierror)
749 stop "Error in double torsional parameter file"
751 read (itordp,*) ntermd_1(i,j,k,iblock),
752 & ntermd_2(i,j,k,iblock)
753 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
754 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
755 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
756 & ntermd_1(i,j,k,iblock))
757 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
758 & ntermd_1(i,j,k,iblock))
759 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
760 & ntermd_1(i,j,k,iblock))
761 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
762 & ntermd_1(i,j,k,iblock))
763 C Martix of D parameters for one dimesional foureir series
764 do l=1,ntermd_1(i,j,k,iblock)
765 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
766 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
767 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
768 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
769 c write(iout,*) "whcodze" ,
770 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
772 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
773 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
774 & v2s(m,l,i,j,k,iblock),
775 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
776 C Martix of D parameters for two dimesional fourier series
777 do l=1,ntermd_2(i,j,k,iblock)
779 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
780 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
781 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
782 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
791 write (iout,*) 'Constants for double torsionals'
794 do j=-ntortyp+1,ntortyp-1
795 do k=-ntortyp+1,ntortyp-1
796 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
797 & ' nsingle',ntermd_1(i,j,k,iblock),
798 & ' ndouble',ntermd_2(i,j,k,iblock)
800 write (iout,*) 'Single angles:'
801 do l=1,ntermd_1(i,j,k,iblock)
802 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
803 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
804 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
805 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
808 write (iout,*) 'Pairs of angles:'
809 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
810 do l=1,ntermd_2(i,j,k,iblock)
811 write (iout,'(i5,20f10.5)')
812 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
815 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
816 do l=1,ntermd_2(i,j,k,iblock)
817 write (iout,'(i5,20f10.5)')
818 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
819 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
828 C Read of Side-chain backbone correlation parameters
829 C Modified 11 May 2012 by Adasko
832 read (isccor,*) nsccortyp
833 read (isccor,*) (isccortyp(i),i=1,ntyp)
835 isccortyp(i)=-isccortyp(-i)
837 iscprol=isccortyp(20)
838 c write (iout,*) 'ntortyp',ntortyp
840 cc maxinter is maximum interaction sites
845 &nterm_sccor(i,j),nlor_sccor(i,j)
846 write (iout,*) nterm_sccor(i,j)
852 nterm_sccor(-i,j)=nterm_sccor(i,j)
853 nterm_sccor(-i,-j)=nterm_sccor(i,j)
854 nterm_sccor(i,-j)=nterm_sccor(i,j)
855 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
856 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
857 do k=1,nterm_sccor(i,j)
858 read (isccor,*) kk,v1sccor(k,l,i,j)
860 if (j.eq.iscprol) then
861 if (i.eq.isccortyp(10)) then
862 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
863 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
865 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
866 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
867 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
868 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
869 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
870 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
871 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
872 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
875 if (i.eq.isccortyp(10)) then
876 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
877 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
879 if (j.eq.isccortyp(10)) then
880 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
881 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
883 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
884 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
885 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
886 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
887 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
888 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
892 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
893 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
894 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
895 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
898 do k=1,nlor_sccor(i,j)
899 read (isccor,*) kk,vlor1sccor(k,i,j),
900 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
901 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
902 &(1+vlor3sccor(k,i,j)**2)
904 v0sccor(l,i,j)=v0ijsccor
905 v0sccor(l,-i,j)=v0ijsccor1
906 v0sccor(l,i,-j)=v0ijsccor2
907 v0sccor(l,-i,-j)=v0ijsccor3
913 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
916 write (iout,*) 'ityp',i,' jtyp',j
917 write (iout,*) 'Fourier constants'
918 do k=1,nterm_sccor(i,j)
919 write (iout,'(2(1pe15.5))')
920 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
922 write (iout,*) 'Lorenz constants'
923 do k=1,nlor_sccor(i,j)
924 write (iout,'(3(1pe15.5))')
925 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
931 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
932 C interaction energy of the Gly, Ala, and Pro prototypes.
934 read (ifourier,*) nloctyp
937 read (ifourier,*) (b(ii,i),ii=1,13)
939 write (iout,*) 'Type',i
940 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
948 B1tilde(1,i) = b(3,i)
949 B1tilde(2,i) =-b(5,i)
950 B1tilde(1,-i) =-b(3,i)
951 B1tilde(2,-i) =b(5,i)
975 Ctilde(2,1,i)=-b(9,i)
977 Ctilde(1,1,-i)=b(7,i)
978 Ctilde(1,2,-i)=-b(9,i)
979 Ctilde(2,1,-i)=b(9,i)
980 Ctilde(2,2,-i)=b(7,i)
982 c Ctilde(1,1,i)=0.0d0
983 c Ctilde(1,2,i)=0.0d0
984 c Ctilde(2,1,i)=0.0d0
985 c Ctilde(2,2,i)=0.0d0
1000 Dtilde(2,1,i)=-b(8,i)
1001 Dtilde(2,2,i)=b(6,i)
1002 Dtilde(1,1,-i)=b(6,i)
1003 Dtilde(1,2,-i)=-b(8,i)
1004 Dtilde(2,1,-i)=b(8,i)
1005 Dtilde(2,2,-i)=b(6,i)
1007 c Dtilde(1,1,i)=0.0d0
1008 c Dtilde(1,2,i)=0.0d0
1009 c Dtilde(2,1,i)=0.0d0
1010 c Dtilde(2,2,i)=0.0d0
1011 EE(1,1,i)= b(10,i)+b(11,i)
1012 EE(2,2,i)=-b(10,i)+b(11,i)
1013 EE(2,1,i)= b(12,i)-b(13,i)
1014 EE(1,2,i)= b(12,i)+b(13,i)
1015 EE(1,1,-i)= b(10,i)+b(11,i)
1016 EE(2,2,-i)=-b(10,i)+b(11,i)
1017 EE(2,1,-i)=-b(12,i)+b(13,i)
1018 EE(1,2,-i)=-b(12,i)-b(13,i)
1024 c ee(2,1,i)=ee(1,2,i)
1029 write (iout,*) 'Type',i
1031 c write (iout,'(f10.5)') B1(:,i)
1032 write(iout,*) B1(1,i),B1(2,i)
1034 c write (iout,'(f10.5)') B2(:,i)
1035 write(iout,*) B2(1,i),B2(2,i)
1038 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1042 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1046 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1051 C Read electrostatic-interaction parameters
1054 write (iout,'(/a)') 'Electrostatic interaction constants:'
1055 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1056 & 'IT','JT','APP','BPP','AEL6','AEL3'
1058 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1059 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1060 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1061 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1066 app (i,j)=epp(i,j)*rri*rri
1067 bpp (i,j)=-2.0D0*epp(i,j)*rri
1068 ael6(i,j)=elpp6(i,j)*4.2D0**6
1069 ael3(i,j)=elpp3(i,j)*4.2D0**3
1071 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1072 & ael6(i,j),ael3(i,j)
1077 C Read side-chain interaction parameters.
1079 read (isidep,*) ipot,expon
1080 if (ipot.lt.1 .or. ipot.gt.5) then
1081 write (iout,'(2a)') 'Error while reading SC interaction',
1082 & 'potential file - unknown potential type.'
1086 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1087 & ', exponents are ',expon,2*expon
1088 goto (10,20,30,30,40) ipot
1089 C----------------------- LJ potential ---------------------------------
1090 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1092 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1093 write (iout,'(a/)') 'The epsilon array:'
1094 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1095 write (iout,'(/a)') 'One-body parameters:'
1096 write (iout,'(a,4x,a)') 'residue','sigma'
1097 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1100 C----------------------- LJK potential --------------------------------
1101 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1102 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1104 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1105 write (iout,'(a/)') 'The epsilon array:'
1106 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1107 write (iout,'(/a)') 'One-body parameters:'
1108 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1109 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1113 C---------------------- GB or BP potential -----------------------------
1115 read (isidep,*)(eps(i,j),j=i,ntyp)
1117 read (isidep,*)(sigma0(i),i=1,ntyp)
1118 read (isidep,*)(sigii(i),i=1,ntyp)
1119 read (isidep,*)(chip(i),i=1,ntyp)
1120 read (isidep,*)(alp(i),i=1,ntyp)
1122 read (isidep,*)(epslip(i,j),j=i,ntyp)
1123 C write(iout,*) "WARNING!!",i,ntyp
1124 write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1126 C epslip(i,j)=epslip(i,j)+0.05d0
1129 C For the GB potential convert sigma'**2 into chi'
1132 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1136 write (iout,'(/a/)') 'Parameters of the BP potential:'
1137 write (iout,'(a/)') 'The epsilon array:'
1138 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1139 write (iout,'(/a)') 'One-body parameters:'
1140 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1142 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1143 & chip(i),alp(i),i=1,ntyp)
1146 C--------------------- GBV potential -----------------------------------
1147 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1148 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1149 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1151 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1152 write (iout,'(a/)') 'The epsilon array:'
1153 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1154 write (iout,'(/a)') 'One-body parameters:'
1155 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1156 & 's||/s_|_^2',' chip ',' alph '
1157 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1158 & sigii(i),chip(i),alp(i),i=1,ntyp)
1162 C-----------------------------------------------------------------------
1163 C Calculate the "working" parameters of SC interactions.
1167 epslip(i,j)=epslip(j,i)
1172 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1173 sigma(j,i)=sigma(i,j)
1174 rs0(i,j)=dwa16*sigma(i,j)
1178 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1179 & 'Working parameters of the SC interactions:',
1180 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1185 epsijlip=epslip(i,j)
1186 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1195 sigeps=dsign(1.0D0,epsij)
1197 aa_aq(i,j)=epsij*rrij*rrij
1198 bb_aq(i,j)=-sigeps*epsij*rrij
1199 aa_aq(j,i)=aa_aq(i,j)
1200 bb_aq(j,i)=bb_aq(i,j)
1201 sigeps=dsign(1.0D0,epsijlip)
1202 epsijlip=dabs(epsijlip)
1203 aa_lip(i,j)=epsijlip*rrij*rrij
1204 bb_lip(i,j)=-sigeps*epsijlip*rrij
1205 aa_lip(j,i)=aa_lip(i,j)
1206 bb_lip(j,i)=bb_lip(i,j)
1208 sigt1sq=sigma0(i)**2
1209 sigt2sq=sigma0(j)**2
1212 ratsig1=sigt2sq/sigt1sq
1213 ratsig2=1.0D0/ratsig1
1214 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1215 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1216 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1220 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1221 sigmaii(i,j)=rsum_max
1222 sigmaii(j,i)=rsum_max
1224 c sigmaii(i,j)=r0(i,j)
1225 c sigmaii(j,i)=r0(i,j)
1227 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1228 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1229 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1230 augm(i,j)=epsij*r_augm**(2*expon)
1231 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1238 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1239 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1240 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1244 write(iout,*) "tube param"
1245 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep,
1246 & ccavtubpep,dcavtubpep,tubetranenepep
1247 sigmapeptube=sigmapeptube**6
1248 sigeps=dsign(1.0D0,epspeptube)
1249 epspeptube=dabs(epspeptube)
1250 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1251 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1252 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1254 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i),
1255 & ccavtub(i),dcavtub(i),tubetranene(i)
1256 sigmasctube=sigmasctube**6
1257 sigeps=dsign(1.0D0,epssctube)
1258 epssctube=dabs(epssctube)
1259 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1260 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1261 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1265 C Define the SC-p interaction constants
1269 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1271 c aad(i,1)=0.3D0*4.0D0**12
1272 C Following line for constants currently implemented
1273 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1274 aad(i,1)=1.5D0*4.0D0**12
1275 c aad(i,1)=0.17D0*5.6D0**12
1277 C "Soft" SC-p repulsion
1279 C Following line for constants currently implemented
1280 c aad(i,1)=0.3D0*4.0D0**6
1281 C "Hard" SC-p repulsion
1282 bad(i,1)=3.0D0*4.0D0**6
1283 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1292 C 8/9/01 Read the SC-p interaction constants from file
1295 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1298 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1299 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1300 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1301 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1305 write (iout,*) "Parameters of SC-p interactions:"
1307 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1308 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1313 C Define the constants of the disulfide bridge
1317 c Old arbitrary potential - commented out.
1322 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1323 c energy surface of diethyl disulfide.
1324 c A. Liwo and U. Kozlowska, 11/24/03
1333 write (iout,*) dyn_ss,'dyndyn'
1335 ss_depth=ebr/wsc-0.25*eps(1,1)
1336 C write(iout,*) akcm,whpb,wsc,'KURWA'
1337 Ht=Ht/wsc-0.25*eps(1,1)
1346 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1350 write (iout,'(/a)') "Disulfide bridge parameters:"
1351 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1352 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1353 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1354 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1357 if (shield_mode.gt.0) then
1359 C VSolvSphere the volume of solving sphere
1361 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1362 C there will be no distinction between proline peptide group and normal peptide
1363 C group in case of shielding parameters
1364 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1365 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1366 write (iout,*) VSolvSphere,VSolvSphere_div
1367 C long axis of side chain
1369 long_r_sidechain(i)=vbldsc0(1,i)
1370 short_r_sidechain(i)=sigma0(i)