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
147 call card_concat(controlcard,.false.)
149 c Return if not own parameters
151 if (iparm.ne.myparm .and. separate_parset) return
153 call reads(controlcard,"BONDPAR",bondname_t,bondname)
154 open (ibond,file=bondname_t,status='old')
156 call reads(controlcard,"THETPAR",thetname_t,thetname)
157 open (ithep,file=thetname_t,status='old')
159 call reads(controlcard,"ROTPAR",rotname_t,rotname)
160 open (irotam,file=rotname_t,status='old')
162 call reads(controlcard,"TORPAR",torname_t,torname)
163 open (itorp,file=torname_t,status='old')
165 call reads(controlcard,"TORDPAR",tordname_t,tordname)
166 open (itordp,file=tordname_t,status='old')
168 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
169 open (isccor,file=sccorname_t,status='old')
171 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
172 open (ifourier,file=fouriername_t,status='old')
174 call reads(controlcard,"ELEPAR",elename_t,elename)
175 open (ielep,file=elename_t,status='old')
177 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
178 open (isidep,file=sidename_t,status='old')
180 call reads(controlcard,"SCPPAR",scpname_t,scpname)
181 open (iscpp,file=scpname_t,status='old')
183 write (iout,*) "Parameter set:",iparm
184 write (iout,*) "Energy-term weights:"
186 write (iout,'(a16,f10.5)') wname(i),ww(i)
188 write (iout,*) "Sidechain potential file : ",
189 & sidename_t(:ilen(sidename_t))
191 write (iout,*) "SCp potential file : ",
192 & scpname_t(:ilen(scpname_t))
194 write (iout,*) "Electrostatic potential file : ",
195 & elename_t(:ilen(elename_t))
196 write (iout,*) "Cumulant coefficient file : ",
197 & fouriername_t(:ilen(fouriername_t))
198 write (iout,*) "Torsional parameter file : ",
199 & torname_t(:ilen(torname_t))
200 write (iout,*) "Double torsional parameter file : ",
201 & tordname_t(:ilen(tordname_t))
202 write (iout,*) "Backbone-rotamer parameter file : ",
203 & sccorname(:ilen(sccorname))
204 write (iout,*) "Bond & inertia constant file : ",
205 & bondname_t(:ilen(bondname_t))
206 write (iout,*) "Bending parameter file : ",
207 & thetname_t(:ilen(thetname_t))
208 write (iout,*) "Rotamer parameter file : ",
209 & rotname_t(:ilen(rotname_t))
212 c Read the virtual-bond parameters, masses, and moments of inertia
213 c and Stokes' radii of the peptide group and side chains
216 read (ibond,*) vbldp0,akp
219 read (ibond,*) vbldsc0(1,i),aksc(1,i)
220 dsc(i) = vbldsc0(1,i)
224 dsc_inv(i)=1.0D0/dsc(i)
228 read (ibond,*) ijunk,vbldp0,akp,rjunk
230 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
232 dsc(i) = vbldsc0(1,i)
236 dsc_inv(i)=1.0D0/dsc(i)
241 write(iout,'(/a/)')"Force constants virtual bonds:"
242 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
244 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
246 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
247 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
249 write (iout,'(13x,3f10.5)')
250 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
256 C Read the parameters of the probability distribution/energy expression
257 C of the virtual-bond valence angles theta
260 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
261 & (bthet(j,i,1,1),j=1,2)
262 read (ithep,*) (polthet(j,i),j=0,3)
263 read (ithep,*) (gthet(j,i),j=1,3)
264 read (ithep,*) theta0(i),sig0(i),sigc0(i)
268 athet(1,i,1,-1)=athet(1,i,1,1)
269 athet(2,i,1,-1)=athet(2,i,1,1)
270 bthet(1,i,1,-1)=-bthet(1,i,1,1)
271 bthet(2,i,1,-1)=-bthet(2,i,1,1)
272 athet(1,i,-1,1)=-athet(1,i,1,1)
273 athet(2,i,-1,1)=-athet(2,i,1,1)
274 bthet(1,i,-1,1)=bthet(1,i,1,1)
275 bthet(2,i,-1,1)=bthet(2,i,1,1)
279 athet(1,i,-1,-1)=athet(1,-i,1,1)
280 athet(2,i,-1,-1)=-athet(2,-i,1,1)
281 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
282 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
283 athet(1,i,-1,1)=athet(1,-i,1,1)
284 athet(2,i,-1,1)=-athet(2,-i,1,1)
285 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
286 bthet(2,i,-1,1)=bthet(2,-i,1,1)
287 athet(1,i,1,-1)=-athet(1,-i,1,1)
288 athet(2,i,1,-1)=athet(2,-i,1,1)
289 bthet(1,i,1,-1)=bthet(1,-i,1,1)
290 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
295 polthet(j,i)=polthet(j,-i)
298 gthet(j,i)=gthet(j,-i)
304 c & 'Parameters of the virtual-bond valence angles:'
305 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
306 c & ' ATHETA0 ',' A1 ',' A2 ',
309 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
310 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
312 c write (iout,'(/a/9x,5a/79(1h-))')
313 c & 'Parameters of the expression for sigma(theta_c):',
314 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
315 c & ' ALPH3 ',' SIGMA0C '
317 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
318 c & (polthet(j,i),j=0,3),sigc0(i)
320 c write (iout,'(/a/9x,5a/79(1h-))')
321 c & 'Parameters of the second gaussian:',
322 c & ' THETA0 ',' SIGMA0 ',' G1 ',
325 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
326 c & sig0(i),(gthet(j,i),j=1,3)
329 & 'Parameters of the virtual-bond valence angles:'
330 write (iout,'(/a/9x,5a/79(1h-))')
331 & 'Coefficients of expansion',
332 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
333 & ' b1*10^1 ',' b2*10^1 '
335 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
336 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
337 & (10*bthet(j,i,1,1),j=1,2)
339 write (iout,'(/a/9x,5a/79(1h-))')
340 & 'Parameters of the expression for sigma(theta_c):',
341 & ' alpha0 ',' alph1 ',' alph2 ',
342 & ' alhp3 ',' sigma0c '
344 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
345 & (polthet(j,i),j=0,3),sigc0(i)
347 write (iout,'(/a/9x,5a/79(1h-))')
348 & 'Parameters of the second gaussian:',
349 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
352 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
353 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
358 C Read the parameters of Utheta determined from ab initio surfaces
359 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
361 write (iout,*) "tu dochodze"
362 read (ithep,*) nthetyp,ntheterm,ntheterm2,
363 & ntheterm3,nsingle,ndouble
364 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
365 read (ithep,*) (ithetyp(i),i=1,ntyp1)
367 ithetyp(i)=-ithetyp(-i)
369 write (iout,*) "tu dochodze"
371 do i=-maxthetyp,maxthetyp
372 do j=-maxthetyp,maxthetyp
373 do k=-maxthetyp,maxthetyp
374 aa0thet(i,j,k,iblock)=0.0d0
376 aathet(l,i,j,k,iblock)=0.0d0
380 bbthet(m,l,i,j,k,iblock)=0.0d0
381 ccthet(m,l,i,j,k,iblock)=0.0d0
382 ddthet(m,l,i,j,k,iblock)=0.0d0
383 eethet(m,l,i,j,k,iblock)=0.0d0
389 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
390 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
398 C write (iout,*) "KURWA1"
401 do j=-nthetyp,nthetyp
402 do k=-nthetyp,nthetyp
403 read (ithep,'(6a)') res1
404 write(iout,*) res1,i,j,k
405 read (ithep,*) aa0thet(i,j,k,iblock)
406 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
408 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
409 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
410 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
411 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
414 & (((ffthet(llll,lll,ll,i,j,k,iblock),
415 & ffthet(lll,llll,ll,i,j,k,iblock),
416 & ggthet(llll,lll,ll,i,j,k,iblock)
417 & ,ggthet(lll,llll,ll,i,j,k,iblock),
418 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
422 C write(iout,*) "KURWA1.1"
424 C For dummy ends assign glycine-type coefficients of theta-only terms; the
425 C coefficients of theta-and-gamma-dependent terms are zero.
430 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
431 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
433 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
434 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
437 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
439 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
442 C write(iout,*) "KURWA1.5"
443 C Substitution for D aminoacids from symmetry.
446 do j=-nthetyp,nthetyp
447 do k=-nthetyp,nthetyp
448 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
450 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
454 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
455 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
456 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
457 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
463 ffthet(llll,lll,ll,i,j,k,iblock)=
464 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
465 ffthet(lll,llll,ll,i,j,k,iblock)=
466 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
467 ggthet(llll,lll,ll,i,j,k,iblock)=
468 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
469 ggthet(lll,llll,ll,i,j,k,iblock)=
470 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
480 C Control printout of the coefficients of virtual-bond-angle potentials
483 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
487 write (iout,'(//4a)')
488 & 'Type ',onelett(i),onelett(j),onelett(k)
489 write (iout,'(//a,10x,a)') " l","a[l]"
490 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
491 write (iout,'(i2,1pe15.5)')
492 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
494 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
495 & "b",l,"c",l,"d",l,"e",l
497 write (iout,'(i2,4(1pe15.5))') m,
498 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
499 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
503 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
504 & "f+",l,"f-",l,"g+",l,"g-",l
507 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
508 & ffthet(n,m,l,i,j,k,iblock),
509 & ffthet(m,n,l,i,j,k,iblock),
510 & ggthet(n,m,l,i,j,k,iblock),
511 & ggthet(m,n,l,i,j,k,iblock)
521 C write(iout,*) 'KURWA2'
524 C Read the parameters of the probability distribution/energy expression
525 C of the side chains.
528 cc write (iout,*) "tu dochodze",i
529 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
533 dsc_inv(i)=1.0D0/dsc(i)
544 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
545 censc(1,1,-i)=censc(1,1,i)
546 censc(2,1,-i)=censc(2,1,i)
547 censc(3,1,-i)=-censc(3,1,i)
549 read (irotam,*) bsc(j,i)
550 read (irotam,*) (censc(k,j,i),k=1,3),
551 & ((blower(k,l,j),l=1,k),k=1,3)
552 censc(1,j,-i)=censc(1,j,i)
553 censc(2,j,-i)=censc(2,j,i)
554 censc(3,j,-i)=-censc(3,j,i)
555 C BSC is amplitude of Gaussian
562 akl=akl+blower(k,m,j)*blower(l,m,j)
566 if (((k.eq.3).and.(l.ne.3))
567 & .or.((l.eq.3).and.(k.ne.3))) then
568 gaussc(k,l,j,-i)=-akl
569 gaussc(l,k,j,-i)=-akl
581 write (iout,'(/a)') 'Parameters of side-chain local geometry'
585 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
586 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
587 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
588 c write (iout,'(a,f10.4,4(16x,f10.4))')
589 c & 'Center ',(bsc(j,i),j=1,nlobi)
590 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
591 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
592 & 'log h',(bsc(j,i),j=1,nlobi)
593 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
594 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
601 c blower(k,l,j)=gaussc(ind,j,i)
606 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
607 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
614 C Read scrot parameters for potentials determined from all-atom AM1 calculations
615 C added by Urszula Kozlowska 07/11/2007
623 read(irotam,*) sc_parmin(j,i)
629 C write (iout,*) 'KURWAKURWA'
632 C Read torsional parameters in old format
634 read (itorp,*) ntortyp,nterm_old
635 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
636 read (itorp,*) (itortyp(i),i=1,ntyp)
641 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
647 write (iout,'(/a/)') 'Torsional constants:'
650 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
651 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
659 C Read torsional parameters
661 read (itorp,*) ntortyp
662 read (itorp,*) (itortyp(i),i=1,ntyp)
663 write (iout,*) 'ntortyp',ntortyp
666 itortyp(i)=-itortyp(-i)
668 c write (iout,*) 'ntortyp',ntortyp
670 do j=-ntortyp+1,ntortyp-1
671 read (itorp,*) nterm(i,j,iblock),
673 nterm(-i,-j,iblock)=nterm(i,j,iblock)
674 nlor(-i,-j,iblock)=nlor(i,j,iblock)
677 do k=1,nterm(i,j,iblock)
678 read (itorp,*) kk,v1(k,i,j,iblock),
680 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
681 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
682 v0ij=v0ij+si*v1(k,i,j,iblock)
685 do k=1,nlor(i,j,iblock)
686 read (itorp,*) kk,vlor1(k,i,j),
687 & vlor2(k,i,j),vlor3(k,i,j)
688 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
691 v0(-i,-j,iblock)=v0ij
697 write (iout,'(/a/)') 'Torsional constants:'
700 write (iout,*) 'ityp',i,' jtyp',j
701 write (iout,*) 'Fourier constants'
702 do k=1,nterm(i,j,iblock)
703 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
706 write (iout,*) 'Lorenz constants'
707 do k=1,nlor(i,j,iblock)
708 write (iout,'(3(1pe15.5))')
709 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
715 C 6/23/01 Read parameters for double torsionals
719 do j=-ntortyp+1,ntortyp-1
720 do k=-ntortyp+1,ntortyp-1
721 read (itordp,'(3a1)') t1,t2,t3
722 c write (iout,*) "OK onelett",
725 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
726 & .or. t3.ne.toronelet(k)) then
727 write (iout,*) "Error in double torsional parameter file",
730 call MPI_Finalize(Ierror)
732 stop "Error in double torsional parameter file"
734 read (itordp,*) ntermd_1(i,j,k,iblock),
735 & ntermd_2(i,j,k,iblock)
736 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
737 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
738 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
739 & ntermd_1(i,j,k,iblock))
740 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
741 & ntermd_1(i,j,k,iblock))
742 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
743 & ntermd_1(i,j,k,iblock))
744 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
745 & ntermd_1(i,j,k,iblock))
746 C Martix of D parameters for one dimesional foureir series
747 do l=1,ntermd_1(i,j,k,iblock)
748 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
749 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
750 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
751 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
752 c write(iout,*) "whcodze" ,
753 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
755 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
756 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
757 & v2s(m,l,i,j,k,iblock),
758 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
759 C Martix of D parameters for two dimesional fourier series
760 do l=1,ntermd_2(i,j,k,iblock)
762 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
763 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
764 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
765 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
774 write (iout,*) 'Constants for double torsionals'
777 do j=-ntortyp+1,ntortyp-1
778 do k=-ntortyp+1,ntortyp-1
779 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
780 & ' nsingle',ntermd_1(i,j,k,iblock),
781 & ' ndouble',ntermd_2(i,j,k,iblock)
783 write (iout,*) 'Single angles:'
784 do l=1,ntermd_1(i,j,k,iblock)
785 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
786 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
787 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
788 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
791 write (iout,*) 'Pairs of angles:'
792 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
793 do l=1,ntermd_2(i,j,k,iblock)
794 write (iout,'(i5,20f10.5)')
795 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
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,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
802 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
811 C Read of Side-chain backbone correlation parameters
812 C Modified 11 May 2012 by Adasko
815 read (isccor,*) nsccortyp
816 read (isccor,*) (isccortyp(i),i=1,ntyp)
818 isccortyp(i)=-isccortyp(-i)
820 iscprol=isccortyp(20)
821 c write (iout,*) 'ntortyp',ntortyp
823 cc maxinter is maximum interaction sites
828 &nterm_sccor(i,j),nlor_sccor(i,j)
829 write (iout,*) nterm_sccor(i,j)
835 nterm_sccor(-i,j)=nterm_sccor(i,j)
836 nterm_sccor(-i,-j)=nterm_sccor(i,j)
837 nterm_sccor(i,-j)=nterm_sccor(i,j)
838 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
839 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
840 do k=1,nterm_sccor(i,j)
841 read (isccor,*) kk,v1sccor(k,l,i,j)
843 if (j.eq.iscprol) then
844 if (i.eq.isccortyp(10)) then
845 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
846 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
848 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
849 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
850 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
851 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
852 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
853 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
854 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
855 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
858 if (i.eq.isccortyp(10)) then
859 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
860 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
862 if (j.eq.isccortyp(10)) then
863 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
864 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
866 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
867 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
868 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
869 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
870 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
871 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
875 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
876 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
877 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
878 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
881 do k=1,nlor_sccor(i,j)
882 read (isccor,*) kk,vlor1sccor(k,i,j),
883 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
884 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
885 &(1+vlor3sccor(k,i,j)**2)
887 v0sccor(l,i,j)=v0ijsccor
888 v0sccor(l,-i,j)=v0ijsccor1
889 v0sccor(l,i,-j)=v0ijsccor2
890 v0sccor(l,-i,-j)=v0ijsccor3
896 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
899 write (iout,*) 'ityp',i,' jtyp',j
900 write (iout,*) 'Fourier constants'
901 do k=1,nterm_sccor(i,j)
902 write (iout,'(2(1pe15.5))')
903 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
905 write (iout,*) 'Lorenz constants'
906 do k=1,nlor_sccor(i,j)
907 write (iout,'(3(1pe15.5))')
908 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
914 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
915 C interaction energy of the Gly, Ala, and Pro prototypes.
917 read (ifourier,*) nloctyp
920 read (ifourier,*) (b(ii,i),ii=1,13)
922 write (iout,*) 'Type',i
923 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
931 B1tilde(1,i) = b(3,i)
932 B1tilde(2,i) =-b(5,i)
933 B1tilde(1,-i) =-b(3,i)
934 B1tilde(2,-i) =b(5,i)
958 Ctilde(2,1,i)=-b(9,i)
960 Ctilde(1,1,-i)=b(7,i)
961 Ctilde(1,2,-i)=-b(9,i)
962 Ctilde(2,1,-i)=b(9,i)
963 Ctilde(2,2,-i)=b(7,i)
965 c Ctilde(1,1,i)=0.0d0
966 c Ctilde(1,2,i)=0.0d0
967 c Ctilde(2,1,i)=0.0d0
968 c Ctilde(2,2,i)=0.0d0
983 Dtilde(2,1,i)=-b(8,i)
985 Dtilde(1,1,-i)=b(6,i)
986 Dtilde(1,2,-i)=-b(8,i)
987 Dtilde(2,1,-i)=b(8,i)
988 Dtilde(2,2,-i)=b(6,i)
990 c Dtilde(1,1,i)=0.0d0
991 c Dtilde(1,2,i)=0.0d0
992 c Dtilde(2,1,i)=0.0d0
993 c Dtilde(2,2,i)=0.0d0
994 EE(1,1,i)= b(10,i)+b(11,i)
995 EE(2,2,i)=-b(10,i)+b(11,i)
996 EE(2,1,i)= b(12,i)-b(13,i)
997 EE(1,2,i)= b(12,i)+b(13,i)
998 EE(1,1,-i)= b(10,i)+b(11,i)
999 EE(2,2,-i)=-b(10,i)+b(11,i)
1000 EE(2,1,-i)=-b(12,i)+b(13,i)
1001 EE(1,2,-i)=-b(12,i)-b(13,i)
1007 c ee(2,1,i)=ee(1,2,i)
1012 write (iout,*) 'Type',i
1014 c write (iout,'(f10.5)') B1(:,i)
1015 write(iout,*) B1(1,i),B1(2,i)
1017 c write (iout,'(f10.5)') B2(:,i)
1018 write(iout,*) B2(1,i),B2(2,i)
1021 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1025 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1029 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1034 C Read electrostatic-interaction parameters
1037 write (iout,'(/a)') 'Electrostatic interaction constants:'
1038 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1039 & 'IT','JT','APP','BPP','AEL6','AEL3'
1041 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1042 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1043 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1044 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1049 app (i,j)=epp(i,j)*rri*rri
1050 bpp (i,j)=-2.0D0*epp(i,j)*rri
1051 ael6(i,j)=elpp6(i,j)*4.2D0**6
1052 ael3(i,j)=elpp3(i,j)*4.2D0**3
1054 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1055 & ael6(i,j),ael3(i,j)
1060 C Read side-chain interaction parameters.
1062 read (isidep,*) ipot,expon
1063 if (ipot.lt.1 .or. ipot.gt.5) then
1064 write (iout,'(2a)') 'Error while reading SC interaction',
1065 & 'potential file - unknown potential type.'
1069 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1070 & ', exponents are ',expon,2*expon
1071 goto (10,20,30,30,40) ipot
1072 C----------------------- LJ potential ---------------------------------
1073 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1075 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1076 write (iout,'(a/)') 'The epsilon array:'
1077 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1078 write (iout,'(/a)') 'One-body parameters:'
1079 write (iout,'(a,4x,a)') 'residue','sigma'
1080 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1083 C----------------------- LJK potential --------------------------------
1084 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1085 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1087 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1088 write (iout,'(a/)') 'The epsilon array:'
1089 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1090 write (iout,'(/a)') 'One-body parameters:'
1091 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1092 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1096 C---------------------- GB or BP potential -----------------------------
1098 read (isidep,*)(eps(i,j),j=i,ntyp)
1100 read (isidep,*)(sigma0(i),i=1,ntyp)
1101 read (isidep,*)(sigii(i),i=1,ntyp)
1102 read (isidep,*)(chip(i),i=1,ntyp)
1103 read (isidep,*)(alp(i),i=1,ntyp)
1104 C For the GB potential convert sigma'**2 into chi'
1107 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1111 write (iout,'(/a/)') 'Parameters of the BP potential:'
1112 write (iout,'(a/)') 'The epsilon array:'
1113 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1114 write (iout,'(/a)') 'One-body parameters:'
1115 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1117 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1118 & chip(i),alp(i),i=1,ntyp)
1121 C--------------------- GBV potential -----------------------------------
1122 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1123 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1124 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1126 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1127 write (iout,'(a/)') 'The epsilon array:'
1128 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1129 write (iout,'(/a)') 'One-body parameters:'
1130 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1131 & 's||/s_|_^2',' chip ',' alph '
1132 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1133 & sigii(i),chip(i),alp(i),i=1,ntyp)
1137 C-----------------------------------------------------------------------
1138 C Calculate the "working" parameters of SC interactions.
1146 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1147 sigma(j,i)=sigma(i,j)
1148 rs0(i,j)=dwa16*sigma(i,j)
1152 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1153 & 'Working parameters of the SC interactions:',
1154 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1159 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1168 sigeps=dsign(1.0D0,epsij)
1170 aa(i,j)=epsij*rrij*rrij
1171 bb(i,j)=-sigeps*epsij*rrij
1175 sigt1sq=sigma0(i)**2
1176 sigt2sq=sigma0(j)**2
1179 ratsig1=sigt2sq/sigt1sq
1180 ratsig2=1.0D0/ratsig1
1181 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1182 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1183 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1187 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1188 sigmaii(i,j)=rsum_max
1189 sigmaii(j,i)=rsum_max
1191 c sigmaii(i,j)=r0(i,j)
1192 c sigmaii(j,i)=r0(i,j)
1194 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1195 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1196 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1197 augm(i,j)=epsij*r_augm**(2*expon)
1198 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1205 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1206 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1207 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1212 C Define the SC-p interaction constants
1216 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1218 c aad(i,1)=0.3D0*4.0D0**12
1219 C Following line for constants currently implemented
1220 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1221 aad(i,1)=1.5D0*4.0D0**12
1222 c aad(i,1)=0.17D0*5.6D0**12
1224 C "Soft" SC-p repulsion
1226 C Following line for constants currently implemented
1227 c aad(i,1)=0.3D0*4.0D0**6
1228 C "Hard" SC-p repulsion
1229 bad(i,1)=3.0D0*4.0D0**6
1230 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1239 C 8/9/01 Read the SC-p interaction constants from file
1242 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1245 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1246 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1247 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1248 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1252 write (iout,*) "Parameters of SC-p interactions:"
1254 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1255 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1260 C Define the constants of the disulfide bridge
1264 c Old arbitrary potential - commented out.
1269 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1270 c energy surface of diethyl disulfide.
1271 c A. Liwo and U. Kozlowska, 11/24/03
1280 write (iout,*) dyn_ss,'dyndyn'
1282 ss_depth=ebr/wsc-0.25*eps(1,1)
1283 C write(iout,*) akcm,whpb,wsc,'KURWA'
1284 Ht=Ht/wsc-0.25*eps(1,1)
1293 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1297 write (iout,'(/a)') "Disulfide bridge parameters:"
1298 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1299 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1300 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1301 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,