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'
25 character*1 onelett(4) /"G","A","P","D"/
26 character*1 toronelet(-2:2) /"p","a","G","A","P"/
28 dimension blower(3,3,maxlob)
29 character*800 controlcard
30 character*256 bondname_t,thetname_t,rotname_t,torname_t,
31 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
37 double precision ip,mp
39 C write (iout,*) "KURWA"
43 C Set LPRINT=.TRUE. for debugging
44 dwa16=2.0d0**(1.0d0/6.0d0)
47 C Assign virtual-bond length
51 call card_concat(controlcard,.true.)
54 key = wname(i)(:ilen(wname(i)))
55 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
58 write (iout,*) "iparm",iparm," myparm",myparm
59 c If reading not own parameters, skip assignment
60 call reada(controlcard,"D0CM",d0cm,3.78d0)
61 call reada(controlcard,"AKCM",akcm,15.1d0)
62 call reada(controlcard,"AKTH",akth,11.0d0)
63 call reada(controlcard,"AKCT",akct,12.0d0)
64 call reada(controlcard,"V1SS",v1ss,-1.08d0)
65 call reada(controlcard,"V2SS",v2ss,7.61d0)
66 call reada(controlcard,"V3SS",v3ss,13.7d0)
67 call reada(controlcard,"EBR",ebr,-5.50D0)
68 call reada(controlcard,"DTRISS",dtriss,1.0D0)
69 call reada(controlcard,"ATRISS",atriss,0.3D0)
70 call reada(controlcard,"BTRISS",btriss,0.02D0)
71 call reada(controlcard,"CTRISS",ctriss,1.0D0)
72 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
73 write(iout,*) "ATRISS",atriss
74 write(iout,*) "BTRISS",btriss
75 write(iout,*) "CTRISS",ctriss
76 write(iout,*) "DTRISS",dtriss
79 C dyn_ss_mask(i)=.false.
83 c Old arbitrary potential - commented out.
88 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
89 c energy surface of diethyl disulfide.
90 c A. Liwo and U. Kozlowska, 11/24/03
102 dyn_ssbond_ij(i,j)=1.0d300
105 call reada(controlcard,"HT",Ht,0.0D0)
107 C ss_depth=ebr/wsc-0.25*eps(1,1)
108 C write(iout,*) HT,wsc,eps(1,1),'KURWA'
109 C Ht=Ht/wsc-0.25*eps(1,1)
118 C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
121 if (iparm.eq.myparm .or. .not.separate_parset) then
124 c Setup weights for UNRES
148 call card_concat(controlcard,.false.)
150 c Return if not own parameters
152 if (iparm.ne.myparm .and. separate_parset) return
154 call reads(controlcard,"BONDPAR",bondname_t,bondname)
155 open (ibond,file=bondname_t,status='old')
157 call reads(controlcard,"THETPAR",thetname_t,thetname)
158 open (ithep,file=thetname_t,status='old')
160 call reads(controlcard,"ROTPAR",rotname_t,rotname)
161 open (irotam,file=rotname_t,status='old')
163 call reads(controlcard,"TORPAR",torname_t,torname)
164 open (itorp,file=torname_t,status='old')
166 call reads(controlcard,"TORDPAR",tordname_t,tordname)
167 open (itordp,file=tordname_t,status='old')
169 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
170 open (isccor,file=sccorname_t,status='old')
172 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
173 open (ifourier,file=fouriername_t,status='old')
175 call reads(controlcard,"ELEPAR",elename_t,elename)
176 open (ielep,file=elename_t,status='old')
178 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
179 open (isidep,file=sidename_t,status='old')
181 call reads(controlcard,"SCPPAR",scpname_t,scpname)
182 open (iscpp,file=scpname_t,status='old')
184 write (iout,*) "Parameter set:",iparm
185 write (iout,*) "Energy-term weights:"
187 write (iout,'(a16,f10.5)') wname(i),ww(i)
189 write (iout,*) "Sidechain potential file : ",
190 & sidename_t(:ilen(sidename_t))
192 write (iout,*) "SCp potential file : ",
193 & scpname_t(:ilen(scpname_t))
195 write (iout,*) "Electrostatic potential file : ",
196 & elename_t(:ilen(elename_t))
197 write (iout,*) "Cumulant coefficient file : ",
198 & fouriername_t(:ilen(fouriername_t))
199 write (iout,*) "Torsional parameter file : ",
200 & torname_t(:ilen(torname_t))
201 write (iout,*) "Double torsional parameter file : ",
202 & tordname_t(:ilen(tordname_t))
203 write (iout,*) "Backbone-rotamer parameter file : ",
204 & sccorname(:ilen(sccorname))
205 write (iout,*) "Bond & inertia constant file : ",
206 & bondname_t(:ilen(bondname_t))
207 write (iout,*) "Bending parameter file : ",
208 & thetname_t(:ilen(thetname_t))
209 write (iout,*) "Rotamer parameter file : ",
210 & rotname_t(:ilen(rotname_t))
213 c Read the virtual-bond parameters, masses, and moments of inertia
214 c and Stokes' radii of the peptide group and side chains
217 read (ibond,*) vbldp0,vbldpdum,akp
220 read (ibond,*) vbldsc0(1,i),aksc(1,i)
221 dsc(i) = vbldsc0(1,i)
225 dsc_inv(i)=1.0D0/dsc(i)
229 read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk
231 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
233 dsc(i) = vbldsc0(1,i)
237 dsc_inv(i)=1.0D0/dsc(i)
242 write(iout,'(/a/)')"Force constants virtual bonds:"
243 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
245 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
247 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
248 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
250 write (iout,'(13x,3f10.5)')
251 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
255 read(iliptranpar,*) pepliptran
257 read(iliptranpar,*) liptranene(i)
262 C Read the parameters of the probability distribution/energy expression
263 C of the virtual-bond valence angles theta
266 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
267 & (bthet(j,i,1,1),j=1,2)
268 read (ithep,*) (polthet(j,i),j=0,3)
269 read (ithep,*) (gthet(j,i),j=1,3)
270 read (ithep,*) theta0(i),sig0(i),sigc0(i)
274 athet(1,i,1,-1)=athet(1,i,1,1)
275 athet(2,i,1,-1)=athet(2,i,1,1)
276 bthet(1,i,1,-1)=-bthet(1,i,1,1)
277 bthet(2,i,1,-1)=-bthet(2,i,1,1)
278 athet(1,i,-1,1)=-athet(1,i,1,1)
279 athet(2,i,-1,1)=-athet(2,i,1,1)
280 bthet(1,i,-1,1)=bthet(1,i,1,1)
281 bthet(2,i,-1,1)=bthet(2,i,1,1)
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)
293 athet(1,i,1,-1)=-athet(1,-i,1,1)
294 athet(2,i,1,-1)=athet(2,-i,1,1)
295 bthet(1,i,1,-1)=bthet(1,-i,1,1)
296 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
301 polthet(j,i)=polthet(j,-i)
304 gthet(j,i)=gthet(j,-i)
310 c & 'Parameters of the virtual-bond valence angles:'
311 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
312 c & ' ATHETA0 ',' A1 ',' A2 ',
315 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
316 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
318 c write (iout,'(/a/9x,5a/79(1h-))')
319 c & 'Parameters of the expression for sigma(theta_c):',
320 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
321 c & ' ALPH3 ',' SIGMA0C '
323 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
324 c & (polthet(j,i),j=0,3),sigc0(i)
326 c write (iout,'(/a/9x,5a/79(1h-))')
327 c & 'Parameters of the second gaussian:',
328 c & ' THETA0 ',' SIGMA0 ',' G1 ',
331 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
332 c & sig0(i),(gthet(j,i),j=1,3)
335 & 'Parameters of the virtual-bond valence angles:'
336 write (iout,'(/a/9x,5a/79(1h-))')
337 & 'Coefficients of expansion',
338 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
339 & ' b1*10^1 ',' b2*10^1 '
341 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
342 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
343 & (10*bthet(j,i,1,1),j=1,2)
345 write (iout,'(/a/9x,5a/79(1h-))')
346 & 'Parameters of the expression for sigma(theta_c):',
347 & ' alpha0 ',' alph1 ',' alph2 ',
348 & ' alhp3 ',' sigma0c '
350 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
351 & (polthet(j,i),j=0,3),sigc0(i)
353 write (iout,'(/a/9x,5a/79(1h-))')
354 & 'Parameters of the second gaussian:',
355 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
358 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
359 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
364 C Read the parameters of Utheta determined from ab initio surfaces
365 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
367 write (iout,*) "tu dochodze"
368 read (ithep,*) nthetyp,ntheterm,ntheterm2,
369 & ntheterm3,nsingle,ndouble
370 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
371 read (ithep,*) (ithetyp(i),i=1,ntyp1)
373 ithetyp(i)=-ithetyp(-i)
375 write (iout,*) "tu dochodze"
377 do i=-maxthetyp,maxthetyp
378 do j=-maxthetyp,maxthetyp
379 do k=-maxthetyp,maxthetyp
380 aa0thet(i,j,k,iblock)=0.0d0
382 aathet(l,i,j,k,iblock)=0.0d0
386 bbthet(m,l,i,j,k,iblock)=0.0d0
387 ccthet(m,l,i,j,k,iblock)=0.0d0
388 ddthet(m,l,i,j,k,iblock)=0.0d0
389 eethet(m,l,i,j,k,iblock)=0.0d0
395 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
396 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
404 C write (iout,*) "KURWA1"
407 do j=-nthetyp,nthetyp
408 do k=-nthetyp,nthetyp
409 read (ithep,'(6a)') res1
410 write(iout,*) res1,i,j,k
411 read (ithep,*) aa0thet(i,j,k,iblock)
412 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
414 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
415 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
416 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
417 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
420 & (((ffthet(llll,lll,ll,i,j,k,iblock),
421 & ffthet(lll,llll,ll,i,j,k,iblock),
422 & ggthet(llll,lll,ll,i,j,k,iblock)
423 & ,ggthet(lll,llll,ll,i,j,k,iblock),
424 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
428 C write(iout,*) "KURWA1.1"
430 C For dummy ends assign glycine-type coefficients of theta-only terms; the
431 C coefficients of theta-and-gamma-dependent terms are zero.
436 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
437 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
439 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
440 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
443 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
445 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
448 C write(iout,*) "KURWA1.5"
449 C Substitution for D aminoacids from symmetry.
452 do j=-nthetyp,nthetyp
453 do k=-nthetyp,nthetyp
454 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
456 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
460 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
461 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
462 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
463 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
469 ffthet(llll,lll,ll,i,j,k,iblock)=
470 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
471 ffthet(lll,llll,ll,i,j,k,iblock)=
472 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
473 ggthet(llll,lll,ll,i,j,k,iblock)=
474 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
475 ggthet(lll,llll,ll,i,j,k,iblock)=
476 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
486 C Control printout of the coefficients of virtual-bond-angle potentials
489 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
493 write (iout,'(//4a)')
494 & 'Type ',onelett(i),onelett(j),onelett(k)
495 write (iout,'(//a,10x,a)') " l","a[l]"
496 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
497 write (iout,'(i2,1pe15.5)')
498 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
500 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
501 & "b",l,"c",l,"d",l,"e",l
503 write (iout,'(i2,4(1pe15.5))') m,
504 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
505 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
509 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
510 & "f+",l,"f-",l,"g+",l,"g-",l
513 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
514 & ffthet(n,m,l,i,j,k,iblock),
515 & ffthet(m,n,l,i,j,k,iblock),
516 & ggthet(n,m,l,i,j,k,iblock),
517 & ggthet(m,n,l,i,j,k,iblock)
527 C write(iout,*) 'KURWA2'
530 C Read the parameters of the probability distribution/energy expression
531 C of the side chains.
534 cc write (iout,*) "tu dochodze",i
535 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
539 dsc_inv(i)=1.0D0/dsc(i)
550 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
551 censc(1,1,-i)=censc(1,1,i)
552 censc(2,1,-i)=censc(2,1,i)
553 censc(3,1,-i)=-censc(3,1,i)
555 read (irotam,*) bsc(j,i)
556 read (irotam,*) (censc(k,j,i),k=1,3),
557 & ((blower(k,l,j),l=1,k),k=1,3)
558 censc(1,j,-i)=censc(1,j,i)
559 censc(2,j,-i)=censc(2,j,i)
560 censc(3,j,-i)=-censc(3,j,i)
561 C BSC is amplitude of Gaussian
568 akl=akl+blower(k,m,j)*blower(l,m,j)
572 if (((k.eq.3).and.(l.ne.3))
573 & .or.((l.eq.3).and.(k.ne.3))) then
574 gaussc(k,l,j,-i)=-akl
575 gaussc(l,k,j,-i)=-akl
587 write (iout,'(/a)') 'Parameters of side-chain local geometry'
591 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
592 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
593 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
594 c write (iout,'(a,f10.4,4(16x,f10.4))')
595 c & 'Center ',(bsc(j,i),j=1,nlobi)
596 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
597 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
598 & 'log h',(bsc(j,i),j=1,nlobi)
599 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
600 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
607 c blower(k,l,j)=gaussc(ind,j,i)
612 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
613 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
620 C Read scrot parameters for potentials determined from all-atom AM1 calculations
621 C added by Urszula Kozlowska 07/11/2007
629 read(irotam,*) sc_parmin(j,i)
635 C write (iout,*) 'KURWAKURWA'
638 C Read torsional parameters in old format
640 read (itorp,*) ntortyp,nterm_old
641 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
642 read (itorp,*) (itortyp(i),i=1,ntyp)
647 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
653 write (iout,'(/a/)') 'Torsional constants:'
656 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
657 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
665 C Read torsional parameters
667 read (itorp,*) ntortyp
668 read (itorp,*) (itortyp(i),i=1,ntyp)
669 write (iout,*) 'ntortyp',ntortyp
672 itortyp(i)=-itortyp(-i)
674 c write (iout,*) 'ntortyp',ntortyp
676 do j=-ntortyp+1,ntortyp-1
677 read (itorp,*) nterm(i,j,iblock),
679 nterm(-i,-j,iblock)=nterm(i,j,iblock)
680 nlor(-i,-j,iblock)=nlor(i,j,iblock)
683 do k=1,nterm(i,j,iblock)
684 read (itorp,*) kk,v1(k,i,j,iblock),
686 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
687 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
688 v0ij=v0ij+si*v1(k,i,j,iblock)
691 do k=1,nlor(i,j,iblock)
692 read (itorp,*) kk,vlor1(k,i,j),
693 & vlor2(k,i,j),vlor3(k,i,j)
694 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
697 v0(-i,-j,iblock)=v0ij
703 write (iout,'(/a/)') 'Torsional constants:'
706 write (iout,*) 'ityp',i,' jtyp',j
707 write (iout,*) 'Fourier constants'
708 do k=1,nterm(i,j,iblock)
709 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
712 write (iout,*) 'Lorenz constants'
713 do k=1,nlor(i,j,iblock)
714 write (iout,'(3(1pe15.5))')
715 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
721 C 6/23/01 Read parameters for double torsionals
725 do j=-ntortyp+1,ntortyp-1
726 do k=-ntortyp+1,ntortyp-1
727 read (itordp,'(3a1)') t1,t2,t3
728 c write (iout,*) "OK onelett",
731 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
732 & .or. t3.ne.toronelet(k)) then
733 write (iout,*) "Error in double torsional parameter file",
736 call MPI_Finalize(Ierror)
738 stop "Error in double torsional parameter file"
740 read (itordp,*) ntermd_1(i,j,k,iblock),
741 & ntermd_2(i,j,k,iblock)
742 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
743 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
744 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
745 & ntermd_1(i,j,k,iblock))
746 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
747 & ntermd_1(i,j,k,iblock))
748 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
749 & ntermd_1(i,j,k,iblock))
750 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
751 & ntermd_1(i,j,k,iblock))
752 C Martix of D parameters for one dimesional foureir series
753 do l=1,ntermd_1(i,j,k,iblock)
754 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
755 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
756 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
757 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
758 c write(iout,*) "whcodze" ,
759 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
761 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
762 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
763 & v2s(m,l,i,j,k,iblock),
764 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
765 C Martix of D parameters for two dimesional fourier series
766 do l=1,ntermd_2(i,j,k,iblock)
768 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
769 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
770 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
771 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
780 write (iout,*) 'Constants for double torsionals'
783 do j=-ntortyp+1,ntortyp-1
784 do k=-ntortyp+1,ntortyp-1
785 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
786 & ' nsingle',ntermd_1(i,j,k,iblock),
787 & ' ndouble',ntermd_2(i,j,k,iblock)
789 write (iout,*) 'Single angles:'
790 do l=1,ntermd_1(i,j,k,iblock)
791 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
792 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
793 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
794 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
797 write (iout,*) 'Pairs of angles:'
798 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
799 do l=1,ntermd_2(i,j,k,iblock)
800 write (iout,'(i5,20f10.5)')
801 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
804 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
805 do l=1,ntermd_2(i,j,k,iblock)
806 write (iout,'(i5,20f10.5)')
807 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
808 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
817 C Read of Side-chain backbone correlation parameters
818 C Modified 11 May 2012 by Adasko
821 read (isccor,*) nsccortyp
822 read (isccor,*) (isccortyp(i),i=1,ntyp)
824 isccortyp(i)=-isccortyp(-i)
826 iscprol=isccortyp(20)
827 c write (iout,*) 'ntortyp',ntortyp
829 cc maxinter is maximum interaction sites
834 &nterm_sccor(i,j),nlor_sccor(i,j)
835 write (iout,*) nterm_sccor(i,j)
841 nterm_sccor(-i,j)=nterm_sccor(i,j)
842 nterm_sccor(-i,-j)=nterm_sccor(i,j)
843 nterm_sccor(i,-j)=nterm_sccor(i,j)
844 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
845 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
846 do k=1,nterm_sccor(i,j)
847 read (isccor,*) kk,v1sccor(k,l,i,j)
849 if (j.eq.iscprol) then
850 if (i.eq.isccortyp(10)) then
851 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
852 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
854 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
855 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
856 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
857 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
858 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
859 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
860 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
861 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
864 if (i.eq.isccortyp(10)) then
865 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
866 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
868 if (j.eq.isccortyp(10)) then
869 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
870 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
872 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
873 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
874 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
875 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
876 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
877 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
881 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
882 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
883 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
884 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
887 do k=1,nlor_sccor(i,j)
888 read (isccor,*) kk,vlor1sccor(k,i,j),
889 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
890 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
891 &(1+vlor3sccor(k,i,j)**2)
893 v0sccor(l,i,j)=v0ijsccor
894 v0sccor(l,-i,j)=v0ijsccor1
895 v0sccor(l,i,-j)=v0ijsccor2
896 v0sccor(l,-i,-j)=v0ijsccor3
902 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
905 write (iout,*) 'ityp',i,' jtyp',j
906 write (iout,*) 'Fourier constants'
907 do k=1,nterm_sccor(i,j)
908 write (iout,'(2(1pe15.5))')
909 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
911 write (iout,*) 'Lorenz constants'
912 do k=1,nlor_sccor(i,j)
913 write (iout,'(3(1pe15.5))')
914 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
920 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
921 C interaction energy of the Gly, Ala, and Pro prototypes.
923 read (ifourier,*) nloctyp
926 read (ifourier,*) (b(ii,i),ii=1,13)
928 write (iout,*) 'Type',i
929 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
937 B1tilde(1,i) = b(3,i)
938 B1tilde(2,i) =-b(5,i)
939 B1tilde(1,-i) =-b(3,i)
940 B1tilde(2,-i) =b(5,i)
964 Ctilde(2,1,i)=-b(9,i)
966 Ctilde(1,1,-i)=b(7,i)
967 Ctilde(1,2,-i)=-b(9,i)
968 Ctilde(2,1,-i)=b(9,i)
969 Ctilde(2,2,-i)=b(7,i)
971 c Ctilde(1,1,i)=0.0d0
972 c Ctilde(1,2,i)=0.0d0
973 c Ctilde(2,1,i)=0.0d0
974 c Ctilde(2,2,i)=0.0d0
989 Dtilde(2,1,i)=-b(8,i)
991 Dtilde(1,1,-i)=b(6,i)
992 Dtilde(1,2,-i)=-b(8,i)
993 Dtilde(2,1,-i)=b(8,i)
994 Dtilde(2,2,-i)=b(6,i)
996 c Dtilde(1,1,i)=0.0d0
997 c Dtilde(1,2,i)=0.0d0
998 c Dtilde(2,1,i)=0.0d0
999 c Dtilde(2,2,i)=0.0d0
1000 EE(1,1,i)= b(10,i)+b(11,i)
1001 EE(2,2,i)=-b(10,i)+b(11,i)
1002 EE(2,1,i)= b(12,i)-b(13,i)
1003 EE(1,2,i)= b(12,i)+b(13,i)
1004 EE(1,1,-i)= b(10,i)+b(11,i)
1005 EE(2,2,-i)=-b(10,i)+b(11,i)
1006 EE(2,1,-i)=-b(12,i)+b(13,i)
1007 EE(1,2,-i)=-b(12,i)-b(13,i)
1013 c ee(2,1,i)=ee(1,2,i)
1018 write (iout,*) 'Type',i
1020 c write (iout,'(f10.5)') B1(:,i)
1021 write(iout,*) B1(1,i),B1(2,i)
1023 c write (iout,'(f10.5)') B2(:,i)
1024 write(iout,*) B2(1,i),B2(2,i)
1027 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1031 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1035 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1040 C Read electrostatic-interaction parameters
1043 write (iout,'(/a)') 'Electrostatic interaction constants:'
1044 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1045 & 'IT','JT','APP','BPP','AEL6','AEL3'
1047 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1048 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1049 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1050 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1055 app (i,j)=epp(i,j)*rri*rri
1056 bpp (i,j)=-2.0D0*epp(i,j)*rri
1057 ael6(i,j)=elpp6(i,j)*4.2D0**6
1058 ael3(i,j)=elpp3(i,j)*4.2D0**3
1060 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1061 & ael6(i,j),ael3(i,j)
1066 C Read side-chain interaction parameters.
1068 read (isidep,*) ipot,expon
1069 if (ipot.lt.1 .or. ipot.gt.5) then
1070 write (iout,'(2a)') 'Error while reading SC interaction',
1071 & 'potential file - unknown potential type.'
1075 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1076 & ', exponents are ',expon,2*expon
1077 goto (10,20,30,30,40) ipot
1078 C----------------------- LJ potential ---------------------------------
1079 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1081 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1082 write (iout,'(a/)') 'The epsilon array:'
1083 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1084 write (iout,'(/a)') 'One-body parameters:'
1085 write (iout,'(a,4x,a)') 'residue','sigma'
1086 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1089 C----------------------- LJK potential --------------------------------
1090 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1091 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1093 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1094 write (iout,'(a/)') 'The epsilon array:'
1095 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1096 write (iout,'(/a)') 'One-body parameters:'
1097 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1098 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1102 C---------------------- GB or BP potential -----------------------------
1104 read (isidep,*)(eps(i,j),j=i,ntyp)
1106 read (isidep,*)(sigma0(i),i=1,ntyp)
1107 read (isidep,*)(sigii(i),i=1,ntyp)
1108 read (isidep,*)(chip(i),i=1,ntyp)
1109 read (isidep,*)(alp(i),i=1,ntyp)
1111 read (isidep,*)(epslip(i,j),j=i,ntyp)
1112 C print *,"WARNING!!"
1114 C epslip(i,j)=epslip(i,j)+0.05d0
1117 C For the GB potential convert sigma'**2 into chi'
1120 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1124 write (iout,'(/a/)') 'Parameters of the BP potential:'
1125 write (iout,'(a/)') 'The epsilon array:'
1126 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1127 write (iout,'(/a)') 'One-body parameters:'
1128 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1130 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1131 & chip(i),alp(i),i=1,ntyp)
1134 C--------------------- GBV potential -----------------------------------
1135 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1136 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1137 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1139 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1140 write (iout,'(a/)') 'The epsilon array:'
1141 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1142 write (iout,'(/a)') 'One-body parameters:'
1143 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1144 & 's||/s_|_^2',' chip ',' alph '
1145 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1146 & sigii(i),chip(i),alp(i),i=1,ntyp)
1150 C-----------------------------------------------------------------------
1151 C Calculate the "working" parameters of SC interactions.
1155 epslip(i,j)=epslip(j,i)
1160 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1161 sigma(j,i)=sigma(i,j)
1162 rs0(i,j)=dwa16*sigma(i,j)
1166 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1167 & 'Working parameters of the SC interactions:',
1168 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1173 epsijlip=epslip(i,j)
1174 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1183 sigeps=dsign(1.0D0,epsij)
1185 aa_aq(i,j)=epsij*rrij*rrij
1186 bb_aq(i,j)=-sigeps*epsij*rrij
1187 aa_aq(j,i)=aa_aq(i,j)
1188 bb_aq(j,i)=bb_aq(i,j)
1189 sigeps=dsign(1.0D0,epsijlip)
1190 epsijlip=dabs(epsijlip)
1191 aa_lip(i,j)=epsijlip*rrij*rrij
1192 bb_lip(i,j)=-sigeps*epsijlip*rrij
1193 aa_lip(j,i)=aa_lip(i,j)
1194 bb_lip(j,i)=bb_lip(i,j)
1196 sigt1sq=sigma0(i)**2
1197 sigt2sq=sigma0(j)**2
1200 ratsig1=sigt2sq/sigt1sq
1201 ratsig2=1.0D0/ratsig1
1202 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1203 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1204 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1208 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1209 sigmaii(i,j)=rsum_max
1210 sigmaii(j,i)=rsum_max
1212 c sigmaii(i,j)=r0(i,j)
1213 c sigmaii(j,i)=r0(i,j)
1215 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1216 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1217 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1218 augm(i,j)=epsij*r_augm**(2*expon)
1219 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1226 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1227 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1228 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1233 C Define the SC-p interaction constants
1237 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1239 c aad(i,1)=0.3D0*4.0D0**12
1240 C Following line for constants currently implemented
1241 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1242 aad(i,1)=1.5D0*4.0D0**12
1243 c aad(i,1)=0.17D0*5.6D0**12
1245 C "Soft" SC-p repulsion
1247 C Following line for constants currently implemented
1248 c aad(i,1)=0.3D0*4.0D0**6
1249 C "Hard" SC-p repulsion
1250 bad(i,1)=3.0D0*4.0D0**6
1251 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1260 C 8/9/01 Read the SC-p interaction constants from file
1263 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1266 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1267 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1268 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1269 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1273 write (iout,*) "Parameters of SC-p interactions:"
1275 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1276 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1281 C Define the constants of the disulfide bridge
1285 c Old arbitrary potential - commented out.
1290 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1291 c energy surface of diethyl disulfide.
1292 c A. Liwo and U. Kozlowska, 11/24/03
1301 write (iout,*) dyn_ss,'dyndyn'
1303 ss_depth=ebr/wsc-0.25*eps(1,1)
1304 C write(iout,*) akcm,whpb,wsc,'KURWA'
1305 Ht=Ht/wsc-0.25*eps(1,1)
1314 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1318 write (iout,'(/a)') "Disulfide bridge parameters:"
1319 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1320 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1321 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1322 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,