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.CONTROL'
26 character*1 onelett(4) /"G","A","P","D"/
27 character*1 toronelet(-2:2) /"p","a","G","A","P"/
29 dimension blower(3,3,maxlob)
30 character*800 controlcard
31 character*256 bondname_t,thetname_t,rotname_t,torname_t,
32 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
38 double precision ip,mp
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
146 call card_concat(controlcard,.false.)
148 c Return if not own parameters
150 if (iparm.ne.myparm .and. separate_parset) return
152 call reads(controlcard,"BONDPAR",bondname_t,bondname)
153 open (ibond,file=bondname_t,status='old')
155 call reads(controlcard,"THETPAR",thetname_t,thetname)
156 open (ithep,file=thetname_t,status='old')
158 call reads(controlcard,"ROTPAR",rotname_t,rotname)
159 open (irotam,file=rotname_t,status='old')
161 call reads(controlcard,"TORPAR",torname_t,torname)
162 open (itorp,file=torname_t,status='old')
164 call reads(controlcard,"TORDPAR",tordname_t,tordname)
165 open (itordp,file=tordname_t,status='old')
167 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
168 open (isccor,file=sccorname_t,status='old')
170 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
171 open (ifourier,file=fouriername_t,status='old')
173 call reads(controlcard,"ELEPAR",elename_t,elename)
174 open (ielep,file=elename_t,status='old')
176 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
177 open (isidep,file=sidename_t,status='old')
179 call reads(controlcard,"SCPPAR",scpname_t,scpname)
180 open (iscpp,file=scpname_t,status='old')
182 write (iout,*) "Parameter set:",iparm
183 write (iout,*) "Energy-term weights:"
185 write (iout,'(a16,f10.5)') wname(i),ww(i)
187 write (iout,*) "Sidechain potential file : ",
188 & sidename_t(:ilen(sidename_t))
190 write (iout,*) "SCp potential file : ",
191 & scpname_t(:ilen(scpname_t))
193 write (iout,*) "Electrostatic potential file : ",
194 & elename_t(:ilen(elename_t))
195 write (iout,*) "Cumulant coefficient file : ",
196 & fouriername_t(:ilen(fouriername_t))
197 write (iout,*) "Torsional parameter file : ",
198 & torname_t(:ilen(torname_t))
199 write (iout,*) "Double torsional parameter file : ",
200 & tordname_t(:ilen(tordname_t))
201 write (iout,*) "Backbone-rotamer parameter file : ",
202 & sccorname(:ilen(sccorname))
203 write (iout,*) "Bond & inertia constant file : ",
204 & bondname_t(:ilen(bondname_t))
205 write (iout,*) "Bending parameter file : ",
206 & thetname_t(:ilen(thetname_t))
207 write (iout,*) "Rotamer parameter file : ",
208 & rotname_t(:ilen(rotname_t))
211 c Read the virtual-bond parameters, masses, and moments of inertia
212 c and Stokes' radii of the peptide group and side chains
215 read (ibond,*) vbldp0,vbldpdum,akp
218 read (ibond,*) vbldsc0(1,i),aksc(1,i)
219 dsc(i) = vbldsc0(1,i)
223 dsc_inv(i)=1.0D0/dsc(i)
227 read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk
229 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
231 dsc(i) = vbldsc0(1,i)
235 dsc_inv(i)=1.0D0/dsc(i)
240 write(iout,'(/a/)')"Force constants virtual bonds:"
241 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
243 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
245 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
246 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
248 write (iout,'(13x,3f10.5)')
249 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
255 C Read the parameters of the probability distribution/energy expression
256 C of the virtual-bond valence angles theta
259 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
260 & (bthet(j,i,1,1),j=1,2)
261 read (ithep,*) (polthet(j,i),j=0,3)
262 read (ithep,*) (gthet(j,i),j=1,3)
263 read (ithep,*) theta0(i),sig0(i),sigc0(i)
267 athet(1,i,1,-1)=athet(1,i,1,1)
268 athet(2,i,1,-1)=athet(2,i,1,1)
269 bthet(1,i,1,-1)=-bthet(1,i,1,1)
270 bthet(2,i,1,-1)=-bthet(2,i,1,1)
271 athet(1,i,-1,1)=-athet(1,i,1,1)
272 athet(2,i,-1,1)=-athet(2,i,1,1)
273 bthet(1,i,-1,1)=bthet(1,i,1,1)
274 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)
282 athet(1,i,-1,1)=athet(1,-i,1,1)
283 athet(2,i,-1,1)=-athet(2,-i,1,1)
284 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
285 bthet(2,i,-1,1)=bthet(2,-i,1,1)
286 athet(1,i,1,-1)=-athet(1,-i,1,1)
287 athet(2,i,1,-1)=athet(2,-i,1,1)
288 bthet(1,i,1,-1)=bthet(1,-i,1,1)
289 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
294 polthet(j,i)=polthet(j,-i)
297 gthet(j,i)=gthet(j,-i)
303 c & 'Parameters of the virtual-bond valence angles:'
304 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
305 c & ' ATHETA0 ',' A1 ',' A2 ',
308 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
309 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
311 c write (iout,'(/a/9x,5a/79(1h-))')
312 c & 'Parameters of the expression for sigma(theta_c):',
313 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
314 c & ' ALPH3 ',' SIGMA0C '
316 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
317 c & (polthet(j,i),j=0,3),sigc0(i)
319 c write (iout,'(/a/9x,5a/79(1h-))')
320 c & 'Parameters of the second gaussian:',
321 c & ' THETA0 ',' SIGMA0 ',' G1 ',
324 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
325 c & sig0(i),(gthet(j,i),j=1,3)
328 & 'Parameters of the virtual-bond valence angles:'
329 write (iout,'(/a/9x,5a/79(1h-))')
330 & 'Coefficients of expansion',
331 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
332 & ' b1*10^1 ',' b2*10^1 '
334 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
335 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
336 & (10*bthet(j,i,1,1),j=1,2)
338 write (iout,'(/a/9x,5a/79(1h-))')
339 & 'Parameters of the expression for sigma(theta_c):',
340 & ' alpha0 ',' alph1 ',' alph2 ',
341 & ' alhp3 ',' sigma0c '
343 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
344 & (polthet(j,i),j=0,3),sigc0(i)
346 write (iout,'(/a/9x,5a/79(1h-))')
347 & 'Parameters of the second gaussian:',
348 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
351 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
352 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
357 C Read the parameters of Utheta determined from ab initio surfaces
358 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
360 read (ithep,*) nthetyp,ntheterm,ntheterm2,
361 & ntheterm3,nsingle,ndouble
362 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
363 read (ithep,*) (ithetyp(i),i=1,ntyp1)
365 ithetyp(i)=-ithetyp(-i)
367 write (iout,*) "tu dochodze"
369 do i=-maxthetyp,maxthetyp
370 do j=-maxthetyp,maxthetyp
371 do k=-maxthetyp,maxthetyp
372 aa0thet(i,j,k,iblock)=0.0d0
374 aathet(l,i,j,k,iblock)=0.0d0
378 bbthet(m,l,i,j,k,iblock)=0.0d0
379 ccthet(m,l,i,j,k,iblock)=0.0d0
380 ddthet(m,l,i,j,k,iblock)=0.0d0
381 eethet(m,l,i,j,k,iblock)=0.0d0
387 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
388 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
396 C write (iout,*) "KURWA1"
399 do j=-nthetyp,nthetyp
400 do k=-nthetyp,nthetyp
401 read (ithep,'(6a)') res1
402 write(iout,*) res1,i,j,k
403 read (ithep,*) aa0thet(i,j,k,iblock)
404 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
406 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
407 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
408 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
409 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
412 & (((ffthet(llll,lll,ll,i,j,k,iblock),
413 & ffthet(lll,llll,ll,i,j,k,iblock),
414 & ggthet(llll,lll,ll,i,j,k,iblock)
415 & ,ggthet(lll,llll,ll,i,j,k,iblock),
416 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
421 C For dummy ends assign glycine-type coefficients of theta-only terms; the
422 C coefficients of theta-and-gamma-dependent terms are zero.
427 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
428 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
430 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
431 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
434 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
436 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
439 C write(iout,*) "KURWA1.5"
440 C Substitution for D aminoacids from symmetry.
443 do j=-nthetyp,nthetyp
444 do k=-nthetyp,nthetyp
445 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
447 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
451 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
452 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
453 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
454 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
460 ffthet(llll,lll,ll,i,j,k,iblock)=
461 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
462 ffthet(lll,llll,ll,i,j,k,iblock)=
463 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
464 ggthet(llll,lll,ll,i,j,k,iblock)=
465 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
466 ggthet(lll,llll,ll,i,j,k,iblock)=
467 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
477 C Control printout of the coefficients of virtual-bond-angle potentials
480 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
484 write (iout,'(//4a)')
485 & 'Type ',onelett(i),onelett(j),onelett(k)
486 write (iout,'(//a,10x,a)') " l","a[l]"
487 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
488 write (iout,'(i2,1pe15.5)')
489 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
491 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
492 & "b",l,"c",l,"d",l,"e",l
494 write (iout,'(i2,4(1pe15.5))') m,
495 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
496 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
500 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
501 & "f+",l,"f-",l,"g+",l,"g-",l
504 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
505 & ffthet(n,m,l,i,j,k,iblock),
506 & ffthet(m,n,l,i,j,k,iblock),
507 & ggthet(n,m,l,i,j,k,iblock),
508 & ggthet(m,n,l,i,j,k,iblock)
521 C Read the parameters of the probability distribution/energy expression
522 C of the side chains.
525 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
529 dsc_inv(i)=1.0D0/dsc(i)
540 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
541 censc(1,1,-i)=censc(1,1,i)
542 censc(2,1,-i)=censc(2,1,i)
543 censc(3,1,-i)=-censc(3,1,i)
545 read (irotam,*) bsc(j,i)
546 read (irotam,*) (censc(k,j,i),k=1,3),
547 & ((blower(k,l,j),l=1,k),k=1,3)
548 censc(1,j,-i)=censc(1,j,i)
549 censc(2,j,-i)=censc(2,j,i)
550 censc(3,j,-i)=-censc(3,j,i)
551 C BSC is amplitude of Gaussian
558 akl=akl+blower(k,m,j)*blower(l,m,j)
562 if (((k.eq.3).and.(l.ne.3))
563 & .or.((l.eq.3).and.(k.ne.3))) then
564 gaussc(k,l,j,-i)=-akl
565 gaussc(l,k,j,-i)=-akl
577 write (iout,'(/a)') 'Parameters of side-chain local geometry'
581 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
582 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
583 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
584 c write (iout,'(a,f10.4,4(16x,f10.4))')
585 c & 'Center ',(bsc(j,i),j=1,nlobi)
586 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
587 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
588 & 'log h',(bsc(j,i),j=1,nlobi)
589 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
590 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
597 c blower(k,l,j)=gaussc(ind,j,i)
602 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
603 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
610 C Read scrot parameters for potentials determined from all-atom AM1 calculations
611 C added by Urszula Kozlowska 07/11/2007
619 read(irotam,*) sc_parmin(j,i)
627 C Read torsional parameters in old format
629 read (itorp,*) ntortyp,nterm_old
630 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
631 read (itorp,*) (itortyp(i),i=1,ntyp)
636 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
642 write (iout,'(/a/)') 'Torsional constants:'
645 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
646 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
654 C Read torsional parameters
656 read (itorp,*) ntortyp
657 read (itorp,*) (itortyp(i),i=1,ntyp)
658 write (iout,*) 'ntortyp',ntortyp
661 itortyp(i)=-itortyp(-i)
663 c write (iout,*) 'ntortyp',ntortyp
665 do j=-ntortyp+1,ntortyp-1
666 read (itorp,*) nterm(i,j,iblock),
668 nterm(-i,-j,iblock)=nterm(i,j,iblock)
669 nlor(-i,-j,iblock)=nlor(i,j,iblock)
672 do k=1,nterm(i,j,iblock)
673 read (itorp,*) kk,v1(k,i,j,iblock),
675 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
676 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
677 v0ij=v0ij+si*v1(k,i,j,iblock)
680 do k=1,nlor(i,j,iblock)
681 read (itorp,*) kk,vlor1(k,i,j),
682 & vlor2(k,i,j),vlor3(k,i,j)
683 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
686 v0(-i,-j,iblock)=v0ij
692 write (iout,'(/a/)') 'Torsional constants:'
695 write (iout,*) 'ityp',i,' jtyp',j
696 write (iout,*) 'Fourier constants'
697 do k=1,nterm(i,j,iblock)
698 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
701 write (iout,*) 'Lorenz constants'
702 do k=1,nlor(i,j,iblock)
703 write (iout,'(3(1pe15.5))')
704 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
710 C 6/23/01 Read parameters for double torsionals
714 do j=-ntortyp+1,ntortyp-1
715 do k=-ntortyp+1,ntortyp-1
716 read (itordp,'(3a1)') t1,t2,t3
717 c write (iout,*) "OK onelett",
720 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
721 & .or. t3.ne.toronelet(k)) then
722 write (iout,*) "Error in double torsional parameter file",
725 call MPI_Finalize(Ierror)
727 stop "Error in double torsional parameter file"
729 read (itordp,*) ntermd_1(i,j,k,iblock),
730 & ntermd_2(i,j,k,iblock)
731 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
732 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
733 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
734 & ntermd_1(i,j,k,iblock))
735 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
736 & ntermd_1(i,j,k,iblock))
737 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
738 & ntermd_1(i,j,k,iblock))
739 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
740 & ntermd_1(i,j,k,iblock))
741 C Martix of D parameters for one dimesional foureir series
742 do l=1,ntermd_1(i,j,k,iblock)
743 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
744 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
745 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
746 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
747 c write(iout,*) "whcodze" ,
748 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
750 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
751 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
752 & v2s(m,l,i,j,k,iblock),
753 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
754 C Martix of D parameters for two dimesional fourier series
755 do l=1,ntermd_2(i,j,k,iblock)
757 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
758 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
759 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
760 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
769 write (iout,*) 'Constants for double torsionals'
772 do j=-ntortyp+1,ntortyp-1
773 do k=-ntortyp+1,ntortyp-1
774 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
775 & ' nsingle',ntermd_1(i,j,k,iblock),
776 & ' ndouble',ntermd_2(i,j,k,iblock)
778 write (iout,*) 'Single angles:'
779 do l=1,ntermd_1(i,j,k,iblock)
780 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
781 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
782 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
783 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
786 write (iout,*) 'Pairs of angles:'
787 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
788 do l=1,ntermd_2(i,j,k,iblock)
789 write (iout,'(i5,20f10.5)')
790 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
793 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
794 do l=1,ntermd_2(i,j,k,iblock)
795 write (iout,'(i5,20f10.5)')
796 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
797 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
806 C Read of Side-chain backbone correlation parameters
807 C Modified 11 May 2012 by Adasko
810 read (isccor,*) nsccortyp
811 read (isccor,*) (isccortyp(i),i=1,ntyp)
813 isccortyp(i)=-isccortyp(-i)
815 iscprol=isccortyp(20)
816 c write (iout,*) 'ntortyp',ntortyp
818 cc maxinter is maximum interaction sites
823 &nterm_sccor(i,j),nlor_sccor(i,j)
824 write (iout,*) nterm_sccor(i,j)
830 nterm_sccor(-i,j)=nterm_sccor(i,j)
831 nterm_sccor(-i,-j)=nterm_sccor(i,j)
832 nterm_sccor(i,-j)=nterm_sccor(i,j)
833 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
834 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
835 do k=1,nterm_sccor(i,j)
836 read (isccor,*) kk,v1sccor(k,l,i,j)
838 if (j.eq.iscprol) then
839 if (i.eq.isccortyp(10)) then
840 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
841 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
843 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
844 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
845 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
846 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
847 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
848 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
849 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
850 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
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 if (j.eq.isccortyp(10)) then
858 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
859 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
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)
865 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
866 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
870 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
871 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
872 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
873 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
876 do k=1,nlor_sccor(i,j)
877 read (isccor,*) kk,vlor1sccor(k,i,j),
878 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
879 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
880 &(1+vlor3sccor(k,i,j)**2)
882 v0sccor(l,i,j)=v0ijsccor
883 v0sccor(l,-i,j)=v0ijsccor1
884 v0sccor(l,i,-j)=v0ijsccor2
885 v0sccor(l,-i,-j)=v0ijsccor3
891 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
894 write (iout,*) 'ityp',i,' jtyp',j
895 write (iout,*) 'Fourier constants'
896 do k=1,nterm_sccor(i,j)
897 write (iout,'(2(1pe15.5))')
898 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
900 write (iout,*) 'Lorenz constants'
901 do k=1,nlor_sccor(i,j)
902 write (iout,'(3(1pe15.5))')
903 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
909 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
910 C interaction energy of the Gly, Ala, and Pro prototypes.
912 read (ifourier,*) nloctyp
915 read (ifourier,*) (b(ii,i),ii=1,13)
917 write (iout,*) 'Type',i
918 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
926 B1tilde(1,i) = b(3,i)
927 B1tilde(2,i) =-b(5,i)
928 B1tilde(1,-i) =-b(3,i)
929 B1tilde(2,-i) =b(5,i)
953 Ctilde(2,1,i)=-b(9,i)
955 Ctilde(1,1,-i)=b(7,i)
956 Ctilde(1,2,-i)=-b(9,i)
957 Ctilde(2,1,-i)=b(9,i)
958 Ctilde(2,2,-i)=b(7,i)
960 c Ctilde(1,1,i)=0.0d0
961 c Ctilde(1,2,i)=0.0d0
962 c Ctilde(2,1,i)=0.0d0
963 c Ctilde(2,2,i)=0.0d0
978 Dtilde(2,1,i)=-b(8,i)
980 Dtilde(1,1,-i)=b(6,i)
981 Dtilde(1,2,-i)=-b(8,i)
982 Dtilde(2,1,-i)=b(8,i)
983 Dtilde(2,2,-i)=b(6,i)
985 c Dtilde(1,1,i)=0.0d0
986 c Dtilde(1,2,i)=0.0d0
987 c Dtilde(2,1,i)=0.0d0
988 c Dtilde(2,2,i)=0.0d0
989 EE(1,1,i)= b(10,i)+b(11,i)
990 EE(2,2,i)=-b(10,i)+b(11,i)
991 EE(2,1,i)= b(12,i)-b(13,i)
992 EE(1,2,i)= b(12,i)+b(13,i)
993 EE(1,1,-i)= b(10,i)+b(11,i)
994 EE(2,2,-i)=-b(10,i)+b(11,i)
995 EE(2,1,-i)=-b(12,i)+b(13,i)
996 EE(1,2,-i)=-b(12,i)-b(13,i)
1002 c ee(2,1,i)=ee(1,2,i)
1007 write (iout,*) 'Type',i
1009 c write (iout,'(f10.5)') B1(:,i)
1010 write(iout,*) B1(1,i),B1(2,i)
1012 c write (iout,'(f10.5)') B2(:,i)
1013 write(iout,*) B2(1,i),B2(2,i)
1016 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1020 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1024 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1029 C Read electrostatic-interaction parameters
1032 write (iout,'(/a)') 'Electrostatic interaction constants:'
1033 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1034 & 'IT','JT','APP','BPP','AEL6','AEL3'
1036 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1037 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1038 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1039 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1044 app (i,j)=epp(i,j)*rri*rri
1045 bpp (i,j)=-2.0D0*epp(i,j)*rri
1046 ael6(i,j)=elpp6(i,j)*4.2D0**6
1047 ael3(i,j)=elpp3(i,j)*4.2D0**3
1049 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1050 & ael6(i,j),ael3(i,j)
1055 C Read side-chain interaction parameters.
1057 read (isidep,*) ipot,expon
1058 if (ipot.lt.1 .or. ipot.gt.5) then
1059 write (iout,'(2a)') 'Error while reading SC interaction',
1060 & 'potential file - unknown potential type.'
1064 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1065 & ', exponents are ',expon,2*expon
1066 goto (10,20,30,30,40) ipot
1067 C----------------------- LJ potential ---------------------------------
1068 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1070 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1071 write (iout,'(a/)') 'The epsilon array:'
1072 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1073 write (iout,'(/a)') 'One-body parameters:'
1074 write (iout,'(a,4x,a)') 'residue','sigma'
1075 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1078 C----------------------- LJK potential --------------------------------
1079 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1080 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1082 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1083 write (iout,'(a/)') 'The epsilon array:'
1084 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1085 write (iout,'(/a)') 'One-body parameters:'
1086 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1087 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1091 C---------------------- GB or BP potential -----------------------------
1093 read (isidep,*)(eps(i,j),j=i,ntyp)
1095 read (isidep,*)(sigma0(i),i=1,ntyp)
1096 read (isidep,*)(sigii(i),i=1,ntyp)
1097 read (isidep,*)(chip(i),i=1,ntyp)
1098 read (isidep,*)(alp(i),i=1,ntyp)
1100 read (isidep,*)(epslip(i,j),j=i,ntyp)
1101 C write(iout,*) "WARNING!!",i,ntyp
1102 write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1104 C epslip(i,j)=epslip(i,j)+0.05d0
1107 C For the GB potential convert sigma'**2 into chi'
1110 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1114 write (iout,'(/a/)') 'Parameters of the BP potential:'
1115 write (iout,'(a/)') 'The epsilon array:'
1116 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1117 write (iout,'(/a)') 'One-body parameters:'
1118 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1120 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1121 & chip(i),alp(i),i=1,ntyp)
1124 C--------------------- GBV potential -----------------------------------
1125 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1126 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1127 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1129 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1130 write (iout,'(a/)') 'The epsilon array:'
1131 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1132 write (iout,'(/a)') 'One-body parameters:'
1133 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1134 & 's||/s_|_^2',' chip ',' alph '
1135 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1136 & sigii(i),chip(i),alp(i),i=1,ntyp)
1140 C-----------------------------------------------------------------------
1141 C Calculate the "working" parameters of SC interactions.
1145 epslip(i,j)=epslip(j,i)
1150 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1151 sigma(j,i)=sigma(i,j)
1152 rs0(i,j)=dwa16*sigma(i,j)
1156 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1157 & 'Working parameters of the SC interactions:',
1158 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1163 epsijlip=epslip(i,j)
1164 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1173 sigeps=dsign(1.0D0,epsij)
1175 aa_aq(i,j)=epsij*rrij*rrij
1176 bb_aq(i,j)=-sigeps*epsij*rrij
1177 aa_aq(j,i)=aa_aq(i,j)
1178 bb_aq(j,i)=bb_aq(i,j)
1179 sigeps=dsign(1.0D0,epsijlip)
1180 epsijlip=dabs(epsijlip)
1181 aa_lip(i,j)=epsijlip*rrij*rrij
1182 bb_lip(i,j)=-sigeps*epsijlip*rrij
1183 aa_lip(j,i)=aa_lip(i,j)
1184 bb_lip(j,i)=bb_lip(i,j)
1186 sigt1sq=sigma0(i)**2
1187 sigt2sq=sigma0(j)**2
1190 ratsig1=sigt2sq/sigt1sq
1191 ratsig2=1.0D0/ratsig1
1192 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1193 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1194 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1198 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1199 sigmaii(i,j)=rsum_max
1200 sigmaii(j,i)=rsum_max
1202 c sigmaii(i,j)=r0(i,j)
1203 c sigmaii(j,i)=r0(i,j)
1205 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1206 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1207 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1208 augm(i,j)=epsij*r_augm**(2*expon)
1209 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1216 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1217 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1218 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1223 C Define the SC-p interaction constants
1227 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1229 c aad(i,1)=0.3D0*4.0D0**12
1230 C Following line for constants currently implemented
1231 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1232 aad(i,1)=1.5D0*4.0D0**12
1233 c aad(i,1)=0.17D0*5.6D0**12
1235 C "Soft" SC-p repulsion
1237 C Following line for constants currently implemented
1238 c aad(i,1)=0.3D0*4.0D0**6
1239 C "Hard" SC-p repulsion
1240 bad(i,1)=3.0D0*4.0D0**6
1241 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1250 C 8/9/01 Read the SC-p interaction constants from file
1253 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1256 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1257 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1258 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1259 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1263 write (iout,*) "Parameters of SC-p interactions:"
1265 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1266 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1271 C Define the constants of the disulfide bridge
1275 c Old arbitrary potential - commented out.
1280 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1281 c energy surface of diethyl disulfide.
1282 c A. Liwo and U. Kozlowska, 11/24/03
1291 write (iout,*) dyn_ss,'dyndyn'
1293 ss_depth=ebr/wsc-0.25*eps(1,1)
1294 C write(iout,*) akcm,whpb,wsc,'KURWA'
1295 Ht=Ht/wsc-0.25*eps(1,1)
1304 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1308 write (iout,'(/a)') "Disulfide bridge parameters:"
1309 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1310 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1311 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1312 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,