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*800 controlcard
32 character*256 bondname_t,thetname_t,rotname_t,torname_t,
33 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_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 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
75 write(iout,*) "ATRISS",atriss
76 write(iout,*) "BTRISS",btriss
77 write(iout,*) "CTRISS",ctriss
78 write(iout,*) "DTRISS",dtriss
81 C dyn_ss_mask(i)=.false.
85 c Old arbitrary potential - commented out.
90 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
91 c energy surface of diethyl disulfide.
92 c A. Liwo and U. Kozlowska, 11/24/03
104 dyn_ssbond_ij(i,j)=1.0d300
107 call reada(controlcard,"HT",Ht,0.0D0)
109 C ss_depth=ebr/wsc-0.25*eps(1,1)
110 C write(iout,*) HT,wsc,eps(1,1),'KURWA'
111 C Ht=Ht/wsc-0.25*eps(1,1)
120 C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
123 if (iparm.eq.myparm .or. .not.separate_parset) then
126 c Setup weights for UNRES
151 call card_concat(controlcard,.false.)
153 c Return if not own parameters
155 if (iparm.ne.myparm .and. separate_parset) return
157 call reads(controlcard,"BONDPAR",bondname_t,bondname)
158 open (ibond,file=bondname_t,status='old')
160 call reads(controlcard,"THETPAR",thetname_t,thetname)
161 open (ithep,file=thetname_t,status='old')
163 call reads(controlcard,"ROTPAR",rotname_t,rotname)
164 open (irotam,file=rotname_t,status='old')
166 call reads(controlcard,"TORPAR",torname_t,torname)
167 open (itorp,file=torname_t,status='old')
169 call reads(controlcard,"TORDPAR",tordname_t,tordname)
170 open (itordp,file=tordname_t,status='old')
172 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
173 open (isccor,file=sccorname_t,status='old')
175 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
176 open (ifourier,file=fouriername_t,status='old')
178 call reads(controlcard,"ELEPAR",elename_t,elename)
179 open (ielep,file=elename_t,status='old')
181 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
182 open (isidep,file=sidename_t,status='old')
184 call reads(controlcard,"SCPPAR",scpname_t,scpname)
185 open (iscpp,file=scpname_t,status='old')
187 write (iout,*) "Parameter set:",iparm
188 write (iout,*) "Energy-term weights:"
190 write (iout,'(a16,f10.5)') wname(i),ww(i)
192 write (iout,*) "Sidechain potential file : ",
193 & sidename_t(:ilen(sidename_t))
195 write (iout,*) "SCp potential file : ",
196 & scpname_t(:ilen(scpname_t))
198 write (iout,*) "Electrostatic potential file : ",
199 & elename_t(:ilen(elename_t))
200 write (iout,*) "Cumulant coefficient file : ",
201 & fouriername_t(:ilen(fouriername_t))
202 write (iout,*) "Torsional parameter file : ",
203 & torname_t(:ilen(torname_t))
204 write (iout,*) "Double torsional parameter file : ",
205 & tordname_t(:ilen(tordname_t))
206 write (iout,*) "Backbone-rotamer parameter file : ",
207 & sccorname(:ilen(sccorname))
208 write (iout,*) "Bond & inertia constant file : ",
209 & bondname_t(:ilen(bondname_t))
210 write (iout,*) "Bending parameter file : ",
211 & thetname_t(:ilen(thetname_t))
212 write (iout,*) "Rotamer parameter file : ",
213 & rotname_t(:ilen(rotname_t))
216 c Read the virtual-bond parameters, masses, and moments of inertia
217 c and Stokes' radii of the peptide group and side chains
220 read (ibond,*) vbldp0,vbldpdum,akp
223 read (ibond,*) vbldsc0(1,i),aksc(1,i)
224 dsc(i) = vbldsc0(1,i)
228 dsc_inv(i)=1.0D0/dsc(i)
232 read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk
234 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
236 dsc(i) = vbldsc0(1,i)
240 dsc_inv(i)=1.0D0/dsc(i)
245 write(iout,'(/a/)')"Force constants virtual bonds:"
246 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
248 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
250 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
251 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
253 write (iout,'(13x,3f10.5)')
254 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
258 read(iliptranpar,*) pepliptran
260 read(iliptranpar,*) liptranene(i)
265 C Read the parameters of the probability distribution/energy expression
266 C of the virtual-bond valence angles theta
269 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
270 & (bthet(j,i,1,1),j=1,2)
271 read (ithep,*) (polthet(j,i),j=0,3)
272 read (ithep,*) (gthet(j,i),j=1,3)
273 read (ithep,*) theta0(i),sig0(i),sigc0(i)
277 athet(1,i,1,-1)=athet(1,i,1,1)
278 athet(2,i,1,-1)=athet(2,i,1,1)
279 bthet(1,i,1,-1)=-bthet(1,i,1,1)
280 bthet(2,i,1,-1)=-bthet(2,i,1,1)
281 athet(1,i,-1,1)=-athet(1,i,1,1)
282 athet(2,i,-1,1)=-athet(2,i,1,1)
283 bthet(1,i,-1,1)=bthet(1,i,1,1)
284 bthet(2,i,-1,1)=bthet(2,i,1,1)
288 athet(1,i,-1,-1)=athet(1,-i,1,1)
289 athet(2,i,-1,-1)=-athet(2,-i,1,1)
290 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
291 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
292 athet(1,i,-1,1)=athet(1,-i,1,1)
293 athet(2,i,-1,1)=-athet(2,-i,1,1)
294 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
295 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)
304 polthet(j,i)=polthet(j,-i)
307 gthet(j,i)=gthet(j,-i)
313 c & 'Parameters of the virtual-bond valence angles:'
314 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
315 c & ' ATHETA0 ',' A1 ',' A2 ',
318 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
319 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
321 c write (iout,'(/a/9x,5a/79(1h-))')
322 c & 'Parameters of the expression for sigma(theta_c):',
323 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
324 c & ' ALPH3 ',' SIGMA0C '
326 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
327 c & (polthet(j,i),j=0,3),sigc0(i)
329 c write (iout,'(/a/9x,5a/79(1h-))')
330 c & 'Parameters of the second gaussian:',
331 c & ' THETA0 ',' SIGMA0 ',' G1 ',
334 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
335 c & sig0(i),(gthet(j,i),j=1,3)
338 & 'Parameters of the virtual-bond valence angles:'
339 write (iout,'(/a/9x,5a/79(1h-))')
340 & 'Coefficients of expansion',
341 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
342 & ' b1*10^1 ',' b2*10^1 '
344 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
345 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
346 & (10*bthet(j,i,1,1),j=1,2)
348 write (iout,'(/a/9x,5a/79(1h-))')
349 & 'Parameters of the expression for sigma(theta_c):',
350 & ' alpha0 ',' alph1 ',' alph2 ',
351 & ' alhp3 ',' sigma0c '
353 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
354 & (polthet(j,i),j=0,3),sigc0(i)
356 write (iout,'(/a/9x,5a/79(1h-))')
357 & 'Parameters of the second gaussian:',
358 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
361 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
362 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
367 C Read the parameters of Utheta determined from ab initio surfaces
368 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
370 write (iout,*) "tu dochodze"
371 read (ithep,*) nthetyp,ntheterm,ntheterm2,
372 & ntheterm3,nsingle,ndouble
373 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
374 read (ithep,*) (ithetyp(i),i=1,ntyp1)
376 ithetyp(i)=-ithetyp(-i)
378 write (iout,*) "tu dochodze"
380 do i=-maxthetyp,maxthetyp
381 do j=-maxthetyp,maxthetyp
382 do k=-maxthetyp,maxthetyp
383 aa0thet(i,j,k,iblock)=0.0d0
385 aathet(l,i,j,k,iblock)=0.0d0
389 bbthet(m,l,i,j,k,iblock)=0.0d0
390 ccthet(m,l,i,j,k,iblock)=0.0d0
391 ddthet(m,l,i,j,k,iblock)=0.0d0
392 eethet(m,l,i,j,k,iblock)=0.0d0
398 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
399 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
407 C write (iout,*) "KURWA1"
410 do j=-nthetyp,nthetyp
411 do k=-nthetyp,nthetyp
412 read (ithep,'(6a)') res1
413 write(iout,*) res1,i,j,k
414 read (ithep,*) aa0thet(i,j,k,iblock)
415 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
417 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
418 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
419 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
420 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
423 & (((ffthet(llll,lll,ll,i,j,k,iblock),
424 & ffthet(lll,llll,ll,i,j,k,iblock),
425 & ggthet(llll,lll,ll,i,j,k,iblock)
426 & ,ggthet(lll,llll,ll,i,j,k,iblock),
427 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
431 C write(iout,*) "KURWA1.1"
433 C For dummy ends assign glycine-type coefficients of theta-only terms; the
434 C coefficients of theta-and-gamma-dependent terms are zero.
439 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
440 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
442 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
443 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
446 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
448 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
451 C write(iout,*) "KURWA1.5"
452 C Substitution for D aminoacids from symmetry.
455 do j=-nthetyp,nthetyp
456 do k=-nthetyp,nthetyp
457 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
459 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
463 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
464 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
465 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
466 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
472 ffthet(llll,lll,ll,i,j,k,iblock)=
473 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
474 ffthet(lll,llll,ll,i,j,k,iblock)=
475 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
476 ggthet(llll,lll,ll,i,j,k,iblock)=
477 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
478 ggthet(lll,llll,ll,i,j,k,iblock)=
479 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
489 C Control printout of the coefficients of virtual-bond-angle potentials
492 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
496 write (iout,'(//4a)')
497 & 'Type ',onelett(i),onelett(j),onelett(k)
498 write (iout,'(//a,10x,a)') " l","a[l]"
499 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
500 write (iout,'(i2,1pe15.5)')
501 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
503 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
504 & "b",l,"c",l,"d",l,"e",l
506 write (iout,'(i2,4(1pe15.5))') m,
507 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
508 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
512 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
513 & "f+",l,"f-",l,"g+",l,"g-",l
516 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
517 & ffthet(n,m,l,i,j,k,iblock),
518 & ffthet(m,n,l,i,j,k,iblock),
519 & ggthet(n,m,l,i,j,k,iblock),
520 & ggthet(m,n,l,i,j,k,iblock)
530 C write(iout,*) 'KURWA2'
533 C Read the parameters of the probability distribution/energy expression
534 C of the side chains.
537 cc write (iout,*) "tu dochodze",i
538 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
542 dsc_inv(i)=1.0D0/dsc(i)
553 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
554 censc(1,1,-i)=censc(1,1,i)
555 censc(2,1,-i)=censc(2,1,i)
556 censc(3,1,-i)=-censc(3,1,i)
558 read (irotam,*) bsc(j,i)
559 read (irotam,*) (censc(k,j,i),k=1,3),
560 & ((blower(k,l,j),l=1,k),k=1,3)
561 censc(1,j,-i)=censc(1,j,i)
562 censc(2,j,-i)=censc(2,j,i)
563 censc(3,j,-i)=-censc(3,j,i)
564 C BSC is amplitude of Gaussian
571 akl=akl+blower(k,m,j)*blower(l,m,j)
575 if (((k.eq.3).and.(l.ne.3))
576 & .or.((l.eq.3).and.(k.ne.3))) then
577 gaussc(k,l,j,-i)=-akl
578 gaussc(l,k,j,-i)=-akl
590 write (iout,'(/a)') 'Parameters of side-chain local geometry'
594 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
595 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
596 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
597 c write (iout,'(a,f10.4,4(16x,f10.4))')
598 c & 'Center ',(bsc(j,i),j=1,nlobi)
599 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
600 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
601 & 'log h',(bsc(j,i),j=1,nlobi)
602 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
603 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
610 c blower(k,l,j)=gaussc(ind,j,i)
615 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
616 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
623 C Read scrot parameters for potentials determined from all-atom AM1 calculations
624 C added by Urszula Kozlowska 07/11/2007
632 read(irotam,*) sc_parmin(j,i)
638 C write (iout,*) 'KURWAKURWA'
641 C Read torsional parameters in old format
643 read (itorp,*) ntortyp,nterm_old
644 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
645 read (itorp,*) (itortyp(i),i=1,ntyp)
650 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
656 write (iout,'(/a/)') 'Torsional constants:'
659 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
660 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
668 C Read torsional parameters
670 read (itorp,*) ntortyp
671 read (itorp,*) (itortyp(i),i=1,ntyp)
672 write (iout,*) 'ntortyp',ntortyp
675 itortyp(i)=-itortyp(-i)
677 c write (iout,*) 'ntortyp',ntortyp
679 do j=-ntortyp+1,ntortyp-1
680 read (itorp,*) nterm(i,j,iblock),
682 nterm(-i,-j,iblock)=nterm(i,j,iblock)
683 nlor(-i,-j,iblock)=nlor(i,j,iblock)
686 do k=1,nterm(i,j,iblock)
687 read (itorp,*) kk,v1(k,i,j,iblock),
689 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
690 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
691 v0ij=v0ij+si*v1(k,i,j,iblock)
694 do k=1,nlor(i,j,iblock)
695 read (itorp,*) kk,vlor1(k,i,j),
696 & vlor2(k,i,j),vlor3(k,i,j)
697 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
700 v0(-i,-j,iblock)=v0ij
706 write (iout,'(/a/)') 'Torsional constants:'
709 write (iout,*) 'ityp',i,' jtyp',j
710 write (iout,*) 'Fourier constants'
711 do k=1,nterm(i,j,iblock)
712 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
715 write (iout,*) 'Lorenz constants'
716 do k=1,nlor(i,j,iblock)
717 write (iout,'(3(1pe15.5))')
718 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
724 C 6/23/01 Read parameters for double torsionals
728 do j=-ntortyp+1,ntortyp-1
729 do k=-ntortyp+1,ntortyp-1
730 read (itordp,'(3a1)') t1,t2,t3
731 c write (iout,*) "OK onelett",
734 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
735 & .or. t3.ne.toronelet(k)) then
736 write (iout,*) "Error in double torsional parameter file",
739 call MPI_Finalize(Ierror)
741 stop "Error in double torsional parameter file"
743 read (itordp,*) ntermd_1(i,j,k,iblock),
744 & ntermd_2(i,j,k,iblock)
745 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
746 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
747 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
748 & ntermd_1(i,j,k,iblock))
749 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
750 & ntermd_1(i,j,k,iblock))
751 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
752 & ntermd_1(i,j,k,iblock))
753 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
754 & ntermd_1(i,j,k,iblock))
755 C Martix of D parameters for one dimesional foureir series
756 do l=1,ntermd_1(i,j,k,iblock)
757 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
758 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
759 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
760 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
761 c write(iout,*) "whcodze" ,
762 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
764 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
765 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
766 & v2s(m,l,i,j,k,iblock),
767 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
768 C Martix of D parameters for two dimesional fourier series
769 do l=1,ntermd_2(i,j,k,iblock)
771 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
772 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
773 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
774 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
783 write (iout,*) 'Constants for double torsionals'
786 do j=-ntortyp+1,ntortyp-1
787 do k=-ntortyp+1,ntortyp-1
788 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
789 & ' nsingle',ntermd_1(i,j,k,iblock),
790 & ' ndouble',ntermd_2(i,j,k,iblock)
792 write (iout,*) 'Single angles:'
793 do l=1,ntermd_1(i,j,k,iblock)
794 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
795 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
796 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
797 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
800 write (iout,*) 'Pairs of angles:'
801 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
802 do l=1,ntermd_2(i,j,k,iblock)
803 write (iout,'(i5,20f10.5)')
804 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
807 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
808 do l=1,ntermd_2(i,j,k,iblock)
809 write (iout,'(i5,20f10.5)')
810 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
811 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
820 C Read of Side-chain backbone correlation parameters
821 C Modified 11 May 2012 by Adasko
824 read (isccor,*) nsccortyp
825 read (isccor,*) (isccortyp(i),i=1,ntyp)
827 isccortyp(i)=-isccortyp(-i)
829 iscprol=isccortyp(20)
830 c write (iout,*) 'ntortyp',ntortyp
832 cc maxinter is maximum interaction sites
837 &nterm_sccor(i,j),nlor_sccor(i,j)
838 write (iout,*) nterm_sccor(i,j)
844 nterm_sccor(-i,j)=nterm_sccor(i,j)
845 nterm_sccor(-i,-j)=nterm_sccor(i,j)
846 nterm_sccor(i,-j)=nterm_sccor(i,j)
847 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
848 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
849 do k=1,nterm_sccor(i,j)
850 read (isccor,*) kk,v1sccor(k,l,i,j)
852 if (j.eq.iscprol) then
853 if (i.eq.isccortyp(10)) then
854 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
855 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
857 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
858 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
859 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
860 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
861 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
862 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
863 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
864 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
867 if (i.eq.isccortyp(10)) then
868 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
869 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
871 if (j.eq.isccortyp(10)) then
872 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
873 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
875 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
876 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
877 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
878 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
879 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
880 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
884 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
885 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
886 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
887 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
890 do k=1,nlor_sccor(i,j)
891 read (isccor,*) kk,vlor1sccor(k,i,j),
892 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
893 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
894 &(1+vlor3sccor(k,i,j)**2)
896 v0sccor(l,i,j)=v0ijsccor
897 v0sccor(l,-i,j)=v0ijsccor1
898 v0sccor(l,i,-j)=v0ijsccor2
899 v0sccor(l,-i,-j)=v0ijsccor3
905 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
908 write (iout,*) 'ityp',i,' jtyp',j
909 write (iout,*) 'Fourier constants'
910 do k=1,nterm_sccor(i,j)
911 write (iout,'(2(1pe15.5))')
912 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
914 write (iout,*) 'Lorenz constants'
915 do k=1,nlor_sccor(i,j)
916 write (iout,'(3(1pe15.5))')
917 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
923 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
924 C interaction energy of the Gly, Ala, and Pro prototypes.
926 read (ifourier,*) nloctyp
929 read (ifourier,*) (b(ii,i),ii=1,13)
931 write (iout,*) 'Type',i
932 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
940 B1tilde(1,i) = b(3,i)
941 B1tilde(2,i) =-b(5,i)
942 B1tilde(1,-i) =-b(3,i)
943 B1tilde(2,-i) =b(5,i)
967 Ctilde(2,1,i)=-b(9,i)
969 Ctilde(1,1,-i)=b(7,i)
970 Ctilde(1,2,-i)=-b(9,i)
971 Ctilde(2,1,-i)=b(9,i)
972 Ctilde(2,2,-i)=b(7,i)
974 c Ctilde(1,1,i)=0.0d0
975 c Ctilde(1,2,i)=0.0d0
976 c Ctilde(2,1,i)=0.0d0
977 c Ctilde(2,2,i)=0.0d0
992 Dtilde(2,1,i)=-b(8,i)
994 Dtilde(1,1,-i)=b(6,i)
995 Dtilde(1,2,-i)=-b(8,i)
996 Dtilde(2,1,-i)=b(8,i)
997 Dtilde(2,2,-i)=b(6,i)
999 c Dtilde(1,1,i)=0.0d0
1000 c Dtilde(1,2,i)=0.0d0
1001 c Dtilde(2,1,i)=0.0d0
1002 c Dtilde(2,2,i)=0.0d0
1003 EE(1,1,i)= b(10,i)+b(11,i)
1004 EE(2,2,i)=-b(10,i)+b(11,i)
1005 EE(2,1,i)= b(12,i)-b(13,i)
1006 EE(1,2,i)= b(12,i)+b(13,i)
1007 EE(1,1,-i)= b(10,i)+b(11,i)
1008 EE(2,2,-i)=-b(10,i)+b(11,i)
1009 EE(2,1,-i)=-b(12,i)+b(13,i)
1010 EE(1,2,-i)=-b(12,i)-b(13,i)
1016 c ee(2,1,i)=ee(1,2,i)
1021 write (iout,*) 'Type',i
1023 c write (iout,'(f10.5)') B1(:,i)
1024 write(iout,*) B1(1,i),B1(2,i)
1026 c write (iout,'(f10.5)') B2(:,i)
1027 write(iout,*) B2(1,i),B2(2,i)
1030 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1034 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1038 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1043 C Read electrostatic-interaction parameters
1046 write (iout,'(/a)') 'Electrostatic interaction constants:'
1047 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1048 & 'IT','JT','APP','BPP','AEL6','AEL3'
1050 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1051 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1052 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1053 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1058 app (i,j)=epp(i,j)*rri*rri
1059 bpp (i,j)=-2.0D0*epp(i,j)*rri
1060 ael6(i,j)=elpp6(i,j)*4.2D0**6
1061 ael3(i,j)=elpp3(i,j)*4.2D0**3
1063 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1064 & ael6(i,j),ael3(i,j)
1069 C Read side-chain interaction parameters.
1071 read (isidep,*) ipot,expon
1072 if (ipot.lt.1 .or. ipot.gt.5) then
1073 write (iout,'(2a)') 'Error while reading SC interaction',
1074 & 'potential file - unknown potential type.'
1078 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1079 & ', exponents are ',expon,2*expon
1080 goto (10,20,30,30,40) ipot
1081 C----------------------- LJ potential ---------------------------------
1082 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1084 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1085 write (iout,'(a/)') 'The epsilon array:'
1086 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1087 write (iout,'(/a)') 'One-body parameters:'
1088 write (iout,'(a,4x,a)') 'residue','sigma'
1089 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1092 C----------------------- LJK potential --------------------------------
1093 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1094 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1096 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1097 write (iout,'(a/)') 'The epsilon array:'
1098 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1099 write (iout,'(/a)') 'One-body parameters:'
1100 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1101 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1105 C---------------------- GB or BP potential -----------------------------
1107 read (isidep,*)(eps(i,j),j=i,ntyp)
1109 read (isidep,*)(sigma0(i),i=1,ntyp)
1110 read (isidep,*)(sigii(i),i=1,ntyp)
1111 read (isidep,*)(chip(i),i=1,ntyp)
1112 read (isidep,*)(alp(i),i=1,ntyp)
1114 read (isidep,*)(epslip(i,j),j=i,ntyp)
1115 C print *,"WARNING!!"
1117 C epslip(i,j)=epslip(i,j)+0.05d0
1120 C For the GB potential convert sigma'**2 into chi'
1123 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1127 write (iout,'(/a/)') 'Parameters of the BP potential:'
1128 write (iout,'(a/)') 'The epsilon array:'
1129 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1130 write (iout,'(/a)') 'One-body parameters:'
1131 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1133 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1134 & chip(i),alp(i),i=1,ntyp)
1137 C--------------------- GBV potential -----------------------------------
1138 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1139 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1140 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1142 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1143 write (iout,'(a/)') 'The epsilon array:'
1144 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1145 write (iout,'(/a)') 'One-body parameters:'
1146 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1147 & 's||/s_|_^2',' chip ',' alph '
1148 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1149 & sigii(i),chip(i),alp(i),i=1,ntyp)
1153 C-----------------------------------------------------------------------
1154 C Calculate the "working" parameters of SC interactions.
1158 epslip(i,j)=epslip(j,i)
1163 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1164 sigma(j,i)=sigma(i,j)
1165 rs0(i,j)=dwa16*sigma(i,j)
1169 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1170 & 'Working parameters of the SC interactions:',
1171 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1176 epsijlip=epslip(i,j)
1177 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1186 sigeps=dsign(1.0D0,epsij)
1188 aa_aq(i,j)=epsij*rrij*rrij
1189 bb_aq(i,j)=-sigeps*epsij*rrij
1190 aa_aq(j,i)=aa_aq(i,j)
1191 bb_aq(j,i)=bb_aq(i,j)
1192 sigeps=dsign(1.0D0,epsijlip)
1193 epsijlip=dabs(epsijlip)
1194 aa_lip(i,j)=epsijlip*rrij*rrij
1195 bb_lip(i,j)=-sigeps*epsijlip*rrij
1196 aa_lip(j,i)=aa_lip(i,j)
1197 bb_lip(j,i)=bb_lip(i,j)
1199 sigt1sq=sigma0(i)**2
1200 sigt2sq=sigma0(j)**2
1203 ratsig1=sigt2sq/sigt1sq
1204 ratsig2=1.0D0/ratsig1
1205 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1206 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1207 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1211 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1212 sigmaii(i,j)=rsum_max
1213 sigmaii(j,i)=rsum_max
1215 c sigmaii(i,j)=r0(i,j)
1216 c sigmaii(j,i)=r0(i,j)
1218 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1219 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1220 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1221 augm(i,j)=epsij*r_augm**(2*expon)
1222 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1229 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1230 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1231 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1236 C Define the SC-p interaction constants
1240 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1242 c aad(i,1)=0.3D0*4.0D0**12
1243 C Following line for constants currently implemented
1244 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1245 aad(i,1)=1.5D0*4.0D0**12
1246 c aad(i,1)=0.17D0*5.6D0**12
1248 C "Soft" SC-p repulsion
1250 C Following line for constants currently implemented
1251 c aad(i,1)=0.3D0*4.0D0**6
1252 C "Hard" SC-p repulsion
1253 bad(i,1)=3.0D0*4.0D0**6
1254 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1263 C 8/9/01 Read the SC-p interaction constants from file
1266 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1269 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1270 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1271 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1272 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1276 write (iout,*) "Parameters of SC-p interactions:"
1278 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1279 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1284 C Define the constants of the disulfide bridge
1288 c Old arbitrary potential - commented out.
1293 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1294 c energy surface of diethyl disulfide.
1295 c A. Liwo and U. Kozlowska, 11/24/03
1304 write (iout,*) dyn_ss,'dyndyn'
1306 ss_depth=ebr/wsc-0.25*eps(1,1)
1307 C write(iout,*) akcm,whpb,wsc,'KURWA'
1308 Ht=Ht/wsc-0.25*eps(1,1)
1317 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1321 write (iout,'(/a)') "Disulfide bridge parameters:"
1322 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1323 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1324 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1325 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1328 if (shield_mode.gt.0) then
1330 C VSolvSphere the volume of solving sphere
1332 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1333 C there will be no distinction between proline peptide group and normal peptide
1334 C group in case of shielding parameters
1335 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1336 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1337 write (iout,*) VSolvSphere,VSolvSphere_div
1338 C long axis of side chain
1340 long_r_sidechain(i)=vbldsc0(1,i)
1341 short_r_sidechain(i)=sigma0(i)