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 C Important! Energy-term weights ARE NOT read here; they are read from the
7 C main input file instead, because NO defaults have yet been set for these
16 include 'COMMON.IOUNITS'
17 include 'COMMON.CHAIN'
18 include 'COMMON.INTERACT'
20 include 'COMMON.LOCAL'
21 include 'COMMON.TORSION'
22 include 'COMMON.SCCOR'
23 include 'COMMON.SCROT'
24 include 'COMMON.FFIELD'
25 include 'COMMON.NAMES'
26 include 'COMMON.SBRIDGE'
29 include 'COMMON.LAGRANGE.5diag'
31 include 'COMMON.LAGRANGE'
35 include 'COMMON.LANGEVIN.lang0.5diag'
37 include 'COMMON.LANGEVIN.lang0'
40 include 'COMMON.LANGEVIN'
42 include 'COMMON.SETUP'
43 include 'COMMON.CONTROL'
44 include 'COMMON.SHIELD'
46 character*1 onelett(4) /"G","A","P","D"/
47 character*1 toronelet(-2:2) /"p","a","G","A","P"/
49 double precision blower(3,3,maxlob)
51 character*3 lancuch,ucase
52 character*1000 weightcard
54 integer i,ii,j,jj,k,kk,l,ll,lll,llll,m,mm,n,iblock,junk,ijunk,
56 double precision akl,v0ij,si,rri,epsij,v0ijsccor,epsijlip,rjunk,
57 & sigt2sq,sigt1sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,
59 double precision dwa16
61 C For printing parameters after they are read set the following in the UNRES
64 C setenv PRINT_PARM YES
66 C To print parameters in LaTeX format rather than as ASCII tables:
70 call getenv_loc("PRINT_PARM",lancuch)
71 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
72 & .and. (me.eq.king.or..not.out1file) .and. fg_rank.eq.0
73 call getenv_loc("LATEX",lancuch)
74 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
76 dwa16=2.0d0**(1.0d0/6.0d0)
79 C Assign virtual-bond length
84 c Read the virtual-bond parameters, masses, and moments of inertia
85 c and Stokes radii of the peptide group and side chains
88 read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
91 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
96 dsc_inv(i)=1.0D0/dsc(i)
100 read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
102 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
103 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
104 dsc(i) = vbldsc0(1,i)
108 dsc_inv(i)=1.0D0/dsc(i)
113 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
114 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
116 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
118 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
119 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
121 write (iout,'(13x,3f10.5)')
122 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
127 C reading lipid parameters
129 write (iout,*) "iliptranpar",iliptranpar
132 read(iliptranpar,*) pepliptran
134 read(iliptranpar,*) liptranene(i)
138 write (iout,'(/a)') "Water-lipid transfer parameters"
139 write (iout,'(a3,3x,f10.5)') 'p',pepliptran
141 write (iout,'(a3,3x,f10.5)') restyp(i),liptranene(i)
146 C Read the parameters of the probability distribution/energy expression
147 C of the virtual-bond valence angles theta
150 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
151 & (bthet(j,i,1,1),j=1,2)
152 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
153 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
154 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
158 athet(1,i,1,-1)=athet(1,i,1,1)
159 athet(2,i,1,-1)=athet(2,i,1,1)
160 bthet(1,i,1,-1)=-bthet(1,i,1,1)
161 bthet(2,i,1,-1)=-bthet(2,i,1,1)
162 athet(1,i,-1,1)=-athet(1,i,1,1)
163 athet(2,i,-1,1)=-athet(2,i,1,1)
164 bthet(1,i,-1,1)=bthet(1,i,1,1)
165 bthet(2,i,-1,1)=bthet(2,i,1,1)
169 athet(1,i,-1,-1)=athet(1,-i,1,1)
170 athet(2,i,-1,-1)=-athet(2,-i,1,1)
171 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
172 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
173 athet(1,i,-1,1)=athet(1,-i,1,1)
174 athet(2,i,-1,1)=-athet(2,-i,1,1)
175 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
176 bthet(2,i,-1,1)=bthet(2,-i,1,1)
177 athet(1,i,1,-1)=-athet(1,-i,1,1)
178 athet(2,i,1,-1)=athet(2,-i,1,1)
179 bthet(1,i,1,-1)=bthet(1,-i,1,1)
180 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
185 polthet(j,i)=polthet(j,-i)
188 gthet(j,i)=gthet(j,-i)
196 & 'Parameters of the virtual-bond valence angles:'
197 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
198 & ' ATHETA0 ',' A1 ',' A2 ',
201 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
202 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
204 write (iout,'(/a/9x,5a/79(1h-))')
205 & 'Parameters of the expression for sigma(theta_c):',
206 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
207 & ' ALPH3 ',' SIGMA0C '
209 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
210 & (polthet(j,i),j=0,3),sigc0(i)
212 write (iout,'(/a/9x,5a/79(1h-))')
213 & 'Parameters of the second gaussian:',
214 & ' THETA0 ',' SIGMA0 ',' G1 ',
217 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
218 & sig0(i),(gthet(j,i),j=1,3)
222 & 'Parameters of the virtual-bond valence angles:'
223 write (iout,'(/a/9x,5a/79(1h-))')
224 & 'Coefficients of expansion',
225 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
226 & ' b1*10^1 ',' b2*10^1 '
228 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
229 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
230 & (10*bthet(j,i,1,1),j=1,2)
232 write (iout,'(/a/9x,5a/79(1h-))')
233 & 'Parameters of the expression for sigma(theta_c):',
234 & ' alpha0 ',' alph1 ',' alph2 ',
235 & ' alhp3 ',' sigma0c '
237 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
238 & (polthet(j,i),j=0,3),sigc0(i)
240 write (iout,'(/a/9x,5a/79(1h-))')
241 & 'Parameters of the second gaussian:',
242 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
245 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
246 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
251 IF (tor_mode.eq.0) THEN
253 C Read the parameters of Utheta determined from ab initio surfaces
254 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
256 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
257 & ntheterm3,nsingle,ndouble
258 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
259 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
261 ithetyp(i)=-ithetyp(-i)
264 do i=-maxthetyp,maxthetyp
265 do j=-maxthetyp,maxthetyp
266 do k=-maxthetyp,maxthetyp
267 aa0thet(i,j,k,iblock)=0.0d0
269 aathet(l,i,j,k,iblock)=0.0d0
273 bbthet(m,l,i,j,k,iblock)=0.0d0
274 ccthet(m,l,i,j,k,iblock)=0.0d0
275 ddthet(m,l,i,j,k,iblock)=0.0d0
276 eethet(m,l,i,j,k,iblock)=0.0d0
282 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
283 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
291 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
293 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
294 c VAR:1=non-glicyne non-proline 2=proline
295 c VAR:negative values for D-aminoacid
297 do j=-nthetyp,nthetyp
298 do k=-nthetyp,nthetyp
299 read (ithep,'(6a)',end=111,err=111) res1
300 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
301 c VAR: aa0thet is variable describing the average value of Foureir
302 c VAR: expansion series
303 c VAR: aathet is foureir expansion in theta/2 angle for full formula
304 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
305 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
306 read (ithep,*,end=111,err=111)
307 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
308 read (ithep,*,end=111,err=111)
309 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
310 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
311 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
312 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
314 read (ithep,*,end=111,err=111)
315 & (((ffthet(llll,lll,ll,i,j,k,iblock),
316 & ffthet(lll,llll,ll,i,j,k,iblock),
317 & ggthet(llll,lll,ll,i,j,k,iblock),
318 & ggthet(lll,llll,ll,i,j,k,iblock),
319 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
324 C For dummy ends assign glycine-type coefficients of theta-only terms; the
325 C coefficients of theta-and-gamma-dependent terms are zero.
326 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
327 C RECOMENTDED AFTER VERSION 3.3)
331 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
332 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
334 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
335 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
338 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
340 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
343 C AND COMMENT THE LOOPS BELOW
347 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
348 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
350 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
351 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
354 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
356 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
360 C Substitution for D aminoacids from symmetry.
363 do j=-nthetyp,nthetyp
364 do k=-nthetyp,nthetyp
365 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
367 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
371 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
372 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
373 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
374 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
380 ffthet(llll,lll,ll,i,j,k,iblock)=
381 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
382 ffthet(lll,llll,ll,i,j,k,iblock)=
383 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
384 ggthet(llll,lll,ll,i,j,k,iblock)=
385 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
386 ggthet(lll,llll,ll,i,j,k,iblock)=
387 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
396 C Control printout of the coefficients of virtual-bond-angle potentials
399 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
402 do j=-nthetyp,nthetyp
403 do k=-nthetyp,nthetyp
404 write (iout,'(//4a)')
405 & 'Type ',toronelet(i),toronelet(j),toronelet(k)
406 write (iout,'(//a,10x,a)') " l","a[l]"
407 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
408 write (iout,'(i2,1pe15.5)')
409 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
411 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
412 & "b",l,"c",l,"d",l,"e",l
414 write (iout,'(i2,4(1pe15.5))') m,
415 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
416 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
420 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
421 & "f+",l,"f-",l,"g+",l,"g-",l
424 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
425 & ffthet(n,m,l,i,j,k,iblock),
426 & ffthet(m,n,l,i,j,k,iblock),
427 & ggthet(n,m,l,i,j,k,iblock),
428 & ggthet(m,n,l,i,j,k,iblock)
441 C here will be the apropriate recalibrating for D-aminoacid
442 read (ithep,*,end=121,err=121) nthetyp
443 do i=-nthetyp+1,nthetyp-1
444 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
445 do j=0,nbend_kcc_Tb(i)
446 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
451 & "Parameters of the valence-only potentials"
452 do i=-nthetyp+1,nthetyp-1
453 write (iout,'(2a)') "Type ",toronelet(i)
454 do k=0,nbend_kcc_Tb(i)
455 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
462 c write (2,*) "Start reading THETA_PDB",ithep_pdb
466 read (ithep_pdb,*,err=111,end=111)
467 & a0thet(i),(athet(j,i,1,1),j=1,2),
468 & (bthet(j,i,1,1),j=1,2)
469 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
470 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
471 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
475 athet(1,i,1,-1)=athet(1,i,1,1)
476 athet(2,i,1,-1)=athet(2,i,1,1)
477 bthet(1,i,1,-1)=-bthet(1,i,1,1)
478 bthet(2,i,1,-1)=-bthet(2,i,1,1)
479 athet(1,i,-1,1)=-athet(1,i,1,1)
480 athet(2,i,-1,1)=-athet(2,i,1,1)
481 bthet(1,i,-1,1)=bthet(1,i,1,1)
482 bthet(2,i,-1,1)=bthet(2,i,1,1)
486 athet(1,i,-1,-1)=athet(1,-i,1,1)
487 athet(2,i,-1,-1)=-athet(2,-i,1,1)
488 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
489 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
490 athet(1,i,-1,1)=athet(1,-i,1,1)
491 athet(2,i,-1,1)=-athet(2,-i,1,1)
492 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
493 bthet(2,i,-1,1)=bthet(2,-i,1,1)
494 athet(1,i,1,-1)=-athet(1,-i,1,1)
495 athet(2,i,1,-1)=athet(2,-i,1,1)
496 bthet(1,i,1,-1)=bthet(1,-i,1,1)
497 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
502 polthet(j,i)=polthet(j,-i)
505 gthet(j,i)=gthet(j,-i)
508 c write (2,*) "End reading THETA_PDB"
514 C Read the parameters of the probability distribution/energy expression
515 C of the side chains.
518 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
522 dsc_inv(i)=1.0D0/dsc(i)
533 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
534 & ((blower(k,l,1),l=1,k),k=1,3)
535 censc(1,1,-i)=censc(1,1,i)
536 censc(2,1,-i)=censc(2,1,i)
537 censc(3,1,-i)=-censc(3,1,i)
539 read (irotam,*,end=112,err=112) bsc(j,i)
540 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
541 & ((blower(k,l,j),l=1,k),k=1,3)
542 censc(1,j,-i)=censc(1,j,i)
543 censc(2,j,-i)=censc(2,j,i)
544 censc(3,j,-i)=-censc(3,j,i)
545 C BSC is amplitude of Gaussian
552 akl=akl+blower(k,m,j)*blower(l,m,j)
556 if (((k.eq.3).and.(l.ne.3))
557 & .or.((l.eq.3).and.(k.ne.3))) then
558 gaussc(k,l,j,-i)=-akl
559 gaussc(l,k,j,-i)=-akl
571 write (iout,'(/a)') 'Parameters of side-chain local geometry'
576 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
577 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
578 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
579 & 'log h',(bsc(j,i),j=1,nlobi)
580 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
581 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
583 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
584 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
587 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
588 write (iout,'(a,f10.4,4(16x,f10.4))')
589 & 'Center ',(bsc(j,i),j=1,nlobi)
590 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
599 C Read scrot parameters for potentials determined from all-atom AM1 calculations
600 C added by Urszula Kozlowska 07/11/2007
603 read (irotam,*,end=112,err=112)
605 read (irotam,*,end=112,err=112)
608 read(irotam,*,end=112,err=112) sc_parmin(j,i)
613 C Read the parameters of the probability distribution/energy expression
614 C of the side chains.
616 write (2,*) "Start reading ROTAM_PDB"
618 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
622 dsc_inv(i)=1.0D0/dsc(i)
633 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
634 & ((blower(k,l,1),l=1,k),k=1,3)
636 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
637 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
638 & ((blower(k,l,j),l=1,k),k=1,3)
645 akl=akl+blower(k,m,j)*blower(l,m,j)
655 write (2,*) "End reading ROTAM_PDB"
659 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
660 C interaction energy of the Gly, Ala, and Pro prototypes.
662 read (ifourier,*) nloctyp
663 SPLIT_FOURIERTOR = nloctyp.lt.0
664 nloctyp = iabs(nloctyp)
666 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
667 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
668 itype2loc(ntyp1)=nloctyp
669 iloctyp(nloctyp)=ntyp1
671 itype2loc(-i)=-itype2loc(i)
680 iloctyp(-i)=-iloctyp(i)
682 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
683 c write (iout,*) "nloctyp",nloctyp,
684 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
693 c write (iout,*) "NEWCORR",i
694 read (ifourier,*,end=115,err=115)
697 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
700 c write (iout,*) "NEWCORR BNEW1"
701 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
704 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
707 c write (iout,*) "NEWCORR BNEW2"
708 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
710 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
711 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
713 c write (iout,*) "NEWCORR CCNEW"
714 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
716 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
717 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
719 c write (iout,*) "NEWCORR DDNEW"
720 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
724 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
728 c write (iout,*) "NEWCORR EENEW1"
729 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
731 read (ifourier,*,end=115,err=115) e0new(ii,i)
733 c write (iout,*) (e0new(ii,i),ii=1,3)
735 c write (iout,*) "NEWCORR EENEW"
738 ccnew(ii,1,i)=ccnew(ii,1,i)/2
739 ccnew(ii,2,i)=ccnew(ii,2,i)/2
740 ddnew(ii,1,i)=ddnew(ii,1,i)/2
741 ddnew(ii,2,i)=ddnew(ii,2,i)/2
746 bnew1(ii,1,-i)= bnew1(ii,1,i)
747 bnew1(ii,2,-i)=-bnew1(ii,2,i)
748 bnew2(ii,1,-i)= bnew2(ii,1,i)
749 bnew2(ii,2,-i)=-bnew2(ii,2,i)
752 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
753 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
754 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
755 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
756 ccnew(ii,1,-i)=ccnew(ii,1,i)
757 ccnew(ii,2,-i)=-ccnew(ii,2,i)
758 ddnew(ii,1,-i)=ddnew(ii,1,i)
759 ddnew(ii,2,-i)=-ddnew(ii,2,i)
761 e0new(1,-i)= e0new(1,i)
762 e0new(2,-i)=-e0new(2,i)
763 e0new(3,-i)=-e0new(3,i)
765 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
766 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
767 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
768 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
772 write (iout,'(a)') "Coefficients of the multibody terms"
773 c do i=-nloctyp+1,nloctyp-1
774 do i=-nloctyp,nloctyp
775 write (iout,*) "Type: ",onelet(iloctyp(i))
776 write (iout,*) "Coefficients of the expansion of B1"
778 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
780 write (iout,*) "Coefficients of the expansion of B2"
782 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
784 write (iout,*) "Coefficients of the expansion of C"
785 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
786 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
787 write (iout,*) "Coefficients of the expansion of D"
788 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
789 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
790 write (iout,*) "Coefficients of the expansion of E"
791 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
794 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
799 IF (SPLIT_FOURIERTOR) THEN
801 c write (iout,*) "NEWCORR TOR",i
802 read (ifourier,*,end=115,err=115)
805 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
808 c write (iout,*) "NEWCORR BNEW1 TOR"
809 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
812 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
815 c write (iout,*) "NEWCORR BNEW2 TOR"
816 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
818 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
819 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
821 c write (iout,*) "NEWCORR CCNEW TOR"
822 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
824 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
825 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
827 c write (iout,*) "NEWCORR DDNEW TOR"
828 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
832 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
836 c write (iout,*) "NEWCORR EENEW1 TOR"
837 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
839 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
841 c write (iout,*) (e0newtor(ii,i),ii=1,3)
843 c write (iout,*) "NEWCORR EENEW TOR"
846 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
847 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
848 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
849 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
854 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
855 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
856 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
857 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
860 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
861 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
862 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
863 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
865 e0newtor(1,-i)= e0newtor(1,i)
866 e0newtor(2,-i)=-e0newtor(2,i)
867 e0newtor(3,-i)=-e0newtor(3,i)
869 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
870 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
871 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
872 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
877 & "Single-body coefficients of the torsional potentials"
878 do i=-nloctyp+1,nloctyp-1
879 write (iout,*) "Type: ",onelet(iloctyp(i))
880 write (iout,*) "Coefficients of the expansion of B1tor"
882 write (iout,'(3hB1(,i1,1h),3f10.5)')
883 & j,(bnew1tor(k,j,i),k=1,3)
885 write (iout,*) "Coefficients of the expansion of B2tor"
887 write (iout,'(3hB2(,i1,1h),3f10.5)')
888 & j,(bnew2tor(k,j,i),k=1,3)
890 write (iout,*) "Coefficients of the expansion of Ctor"
891 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
892 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
893 write (iout,*) "Coefficients of the expansion of Dtor"
894 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
895 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
896 write (iout,*) "Coefficients of the expansion of Etor"
897 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
900 write (iout,'(1hE,2i1,2f10.5)')
901 & j,k,(eenewtor(l,j,k,i),l=1,2)
907 do i=-nloctyp+1,nloctyp-1
910 bnew1tor(ii,j,i)=bnew1(ii,j,i)
915 bnew2tor(ii,j,i)=bnew2(ii,j,i)
919 ccnewtor(ii,1,i)=ccnew(ii,1,i)
920 ccnewtor(ii,2,i)=ccnew(ii,2,i)
921 ddnewtor(ii,1,i)=ddnew(ii,1,i)
922 ddnewtor(ii,2,i)=ddnew(ii,2,i)
925 ENDIF !SPLIT_FOURIER_TOR
928 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
930 read (ifourier,*,end=115,err=115)
931 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
933 write (iout,*) 'Type ',onelet(iloctyp(i))
934 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
948 c B1tilde(1,i) = b(3)
949 c B1tilde(2,i) =-b(5)
950 c B1tilde(1,-i) =-b(3)
951 c B1tilde(2,-i) =b(5)
958 cc B1tilde(1,i) = b(3,i)
959 cc B1tilde(2,i) =-b(5,i)
960 C B1tilde(1,-i) =-b(3,i)
961 C B1tilde(2,-i) =b(5,i)
962 cc b1tilde(1,i)=0.0d0
963 cc b1tilde(2,i)=0.0d0
975 CCold(1,1,-i)= b(7,i)
976 CCold(2,2,-i)=-b(7,i)
977 CCold(2,1,-i)=-b(9,i)
978 CCold(1,2,-i)=-b(9,i)
983 c Ctilde(1,1,i)= CCold(1,1,i)
984 c Ctilde(1,2,i)= CCold(1,2,i)
985 c Ctilde(2,1,i)=-CCold(2,1,i)
986 c Ctilde(2,2,i)=-CCold(2,2,i)
988 c Ctilde(1,1,i)=0.0d0
989 c Ctilde(1,2,i)=0.0d0
990 c Ctilde(2,1,i)=0.0d0
991 c Ctilde(2,2,i)=0.0d0
996 DDold(1,1,-i)= b(6,i)
997 DDold(2,2,-i)=-b(6,i)
998 DDold(2,1,-i)=-b(8,i)
999 DDold(1,2,-i)=-b(8,i)
1004 c Dtilde(1,1,i)= DD(1,1,i)
1005 c Dtilde(1,2,i)= DD(1,2,i)
1006 c Dtilde(2,1,i)=-DD(2,1,i)
1007 c Dtilde(2,2,i)=-DD(2,2,i)
1009 c Dtilde(1,1,i)=0.0d0
1010 c Dtilde(1,2,i)=0.0d0
1011 c Dtilde(2,1,i)=0.0d0
1012 c Dtilde(2,2,i)=0.0d0
1013 EEold(1,1,i)= b(10,i)+b(11,i)
1014 EEold(2,2,i)=-b(10,i)+b(11,i)
1015 EEold(2,1,i)= b(12,i)-b(13,i)
1016 EEold(1,2,i)= b(12,i)+b(13,i)
1017 EEold(1,1,-i)= b(10,i)+b(11,i)
1018 EEold(2,2,-i)=-b(10,i)+b(11,i)
1019 EEold(2,1,-i)=-b(12,i)+b(13,i)
1020 EEold(1,2,-i)=-b(12,i)-b(13,i)
1021 write(iout,*) "TU DOCHODZE"
1027 c ee(2,1,i)=ee(1,2,i)
1032 &"Coefficients of the cumulants (independent of valence angles)"
1033 do i=-nloctyp+1,nloctyp-1
1034 write (iout,*) 'Type ',onelet(iloctyp(i))
1036 write(iout,'(2f10.5)') B(3,i),B(5,i)
1038 write(iout,'(2f10.5)') B(2,i),B(4,i)
1041 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1045 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1049 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1057 C Read torsional parameters in old format
1059 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1060 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1061 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1066 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1072 write (iout,'(/a/)') 'Torsional constants:'
1075 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1076 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1082 C Read torsional parameters
1084 IF (TOR_MODE.eq.0) THEN
1086 read (itorp,*,end=113,err=113) ntortyp
1087 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1089 itype2loc(i)=itortyp(i)
1092 itype2loc(-i)=-itype2loc(i)
1094 itortyp(ntyp1)=ntortyp
1097 itortyp(i)=-itortyp(-i)
1099 write (iout,*) 'ntortyp',ntortyp
1101 do j=-ntortyp+1,ntortyp-1
1102 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1104 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1105 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1108 do k=1,nterm(i,j,iblock)
1109 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1111 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1112 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1113 v0ij=v0ij+si*v1(k,i,j,iblock)
1115 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
1116 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
1117 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
1119 do k=1,nlor(i,j,iblock)
1120 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1121 & vlor2(k,i,j),vlor3(k,i,j)
1122 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1125 v0(-i,-j,iblock)=v0ij
1131 write (iout,'(/a/)') 'Parameters of the SCCOR potentials:'
1134 do j=-ntortyp+1,ntortyp-1
1135 write (iout,*) 'ityp',i,' jtyp',j
1136 write (iout,*) 'Fourier constants'
1137 do k=1,nterm(i,j,iblock)
1138 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1141 write (iout,*) 'Lorenz constants'
1142 do k=1,nlor(i,j,iblock)
1143 write (iout,'(3(1pe15.5))')
1144 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1152 C 6/23/01 Read parameters for double torsionals
1156 do j=-ntortyp+1,ntortyp-1
1157 do k=-ntortyp+1,ntortyp-1
1158 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1159 c write (iout,*) "OK onelett",
1162 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1163 & .or. t3.ne.toronelet(k)) then
1164 write (iout,*) "Error in double torsional parameter file",
1167 call MPI_Finalize(Ierror)
1169 stop "Error in double torsional parameter file"
1171 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1172 & ntermd_2(i,j,k,iblock)
1173 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1174 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1175 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1176 & ntermd_1(i,j,k,iblock))
1177 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1178 & ntermd_1(i,j,k,iblock))
1179 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1180 & ntermd_1(i,j,k,iblock))
1181 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1182 & ntermd_1(i,j,k,iblock))
1183 C Martix of D parameters for one dimesional foureir series
1184 do l=1,ntermd_1(i,j,k,iblock)
1185 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1186 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1187 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1188 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1189 c write(iout,*) "whcodze" ,
1190 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1192 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1193 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1194 & v2s(m,l,i,j,k,iblock),
1195 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1196 C Martix of D parameters for two dimesional fourier series
1197 do l=1,ntermd_2(i,j,k,iblock)
1199 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1200 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1201 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1202 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1211 write (iout,*) 'Constants for double torsionals'
1214 do j=-ntortyp+1,ntortyp-1
1215 do k=-ntortyp+1,ntortyp-1
1216 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1217 & ' nsingle',ntermd_1(i,j,k,iblock),
1218 & ' ndouble',ntermd_2(i,j,k,iblock)
1220 write (iout,*) 'Single angles:'
1221 do l=1,ntermd_1(i,j,k,iblock)
1222 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1223 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1224 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1225 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1228 write (iout,*) 'Pairs of angles:'
1229 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1230 do l=1,ntermd_2(i,j,k,iblock)
1231 write (iout,'(i5,20f10.5)')
1232 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1235 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1236 do l=1,ntermd_2(i,j,k,iblock)
1237 write (iout,'(i5,20f10.5)')
1238 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1239 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1248 ELSE IF (TOR_MODE.eq.1) THEN
1250 C read valence-torsional parameters
1251 read (itorp,*,end=121,err=121) ntortyp
1253 write (iout,*) "Valence-torsional parameters read in ntortyp",
1255 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
1256 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1259 itype2loc(i)=itortyp(i)
1263 itortyp(i)=-itortyp(-i)
1265 do i=-ntortyp+1,ntortyp-1
1266 do j=-ntortyp+1,ntortyp-1
1267 C first we read the cos and sin gamma parameters
1268 read (itorp,'(13x,a)',end=121,err=121) string
1269 write (iout,*) i,j,string
1270 read (itorp,*,end=121,err=121)
1271 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1272 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1273 do k=1,nterm_kcc(j,i)
1274 do l=1,nterm_kcc_Tb(j,i)
1275 do ll=1,nterm_kcc_Tb(j,i)
1276 read (itorp,*,end=121,err=121) ii,jj,kk,
1277 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1285 c AL 4/8/16: Calculate coefficients from one-body parameters
1288 itortyp(i)=itype2loc(i)
1291 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1293 do i=-ntortyp+1,ntortyp-1
1294 do j=-ntortyp+1,ntortyp-1
1297 do k=1,nterm_kcc_Tb(j,i)
1298 do l=1,nterm_kcc_Tb(j,i)
1299 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1300 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1301 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1302 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1305 do k=1,nterm_kcc_Tb(j,i)
1306 do l=1,nterm_kcc_Tb(j,i)
1308 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1309 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1310 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1311 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1313 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1314 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1315 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1316 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1320 c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
1324 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1329 if (tor_mode.gt.0 .and. lprint) then
1330 c Print valence-torsional parameters
1332 & "Parameters of the valence-torsional potentials"
1333 do i=-ntortyp+1,ntortyp-1
1334 do j=-ntortyp+1,ntortyp-1
1335 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1336 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1337 do k=1,nterm_kcc(j,i)
1338 do l=1,nterm_kcc_Tb(j,i)
1339 do ll=1,nterm_kcc_Tb(j,i)
1340 write (iout,'(3i5,2f15.4)')
1341 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1350 C Read of Side-chain backbone correlation parameters
1351 C Modified 11 May 2012 by Adasko
1354 read (isccor,*,end=119,err=119) nsccortyp
1356 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1358 isccortyp(i)=-isccortyp(-i)
1360 iscprol=isccortyp(20)
1361 c write (iout,*) 'ntortyp',ntortyp
1363 cc maxinter is maximum interaction sites
1367 read (isccor,*,end=119,err=119)
1368 &nterm_sccor(i,j),nlor_sccor(i,j)
1374 nterm_sccor(-i,j)=nterm_sccor(i,j)
1375 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1376 nterm_sccor(i,-j)=nterm_sccor(i,j)
1377 do k=1,nterm_sccor(i,j)
1378 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1380 if (j.eq.iscprol) then
1381 if (i.eq.isccortyp(10)) then
1382 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1383 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1385 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1386 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1387 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1388 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1389 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1390 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1391 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1392 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1395 if (i.eq.isccortyp(10)) then
1396 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1397 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1399 if (j.eq.isccortyp(10)) then
1400 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1401 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1403 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1404 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1405 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1406 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1407 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1408 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1412 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1413 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1414 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1415 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1418 do k=1,nlor_sccor(i,j)
1419 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1420 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1421 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1422 &(1+vlor3sccor(k,i,j)**2)
1424 v0sccor(l,i,j)=v0ijsccor
1425 v0sccor(l,-i,j)=v0ijsccor1
1426 v0sccor(l,i,-j)=v0ijsccor2
1427 v0sccor(l,-i,-j)=v0ijsccor3
1433 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1434 c write (iout,*) 'ntortyp',ntortyp
1436 cc maxinter is maximum interaction sites
1440 read (isccor,*,end=119,err=119)
1441 & nterm_sccor(i,j),nlor_sccor(i,j)
1445 do k=1,nterm_sccor(i,j)
1446 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1448 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1451 do k=1,nlor_sccor(i,j)
1452 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1453 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1454 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1455 &(1+vlor3sccor(k,i,j)**2)
1457 v0sccor(l,i,j)=v0ijsccor
1465 write (iout,'(/a/)') 'SCCor torsional constants:'
1469 write (iout,*) 'ityp',i,' jtyp',j
1470 write (iout,*) 'Fourier constants'
1471 do k=1,nterm_sccor(i,j)
1472 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1474 write (iout,*) 'Lorenz constants'
1475 do k=1,nlor_sccor(i,j)
1476 write (iout,'(3(1pe15.5))')
1477 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1486 C Read electrostatic-interaction parameters
1490 write (iout,'(/a)') 'Electrostatic interaction constants:'
1491 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1492 & 'IT','JT','APP','BPP','AEL6','AEL3'
1494 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1495 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1496 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1497 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1502 app (i,j)=epp(i,j)*rri*rri
1503 bpp (i,j)=-2.0D0*epp(i,j)*rri
1504 ael6(i,j)=elpp6(i,j)*4.2D0**6
1505 ael3(i,j)=elpp3(i,j)*4.2D0**3
1507 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1508 & ael6(i,j),ael3(i,j)
1513 C Read side-chain interaction parameters.
1515 read (isidep,*,end=117,err=117) ipot,expon
1516 if (ipot.lt.1 .or. ipot.gt.5) then
1517 write (iout,'(2a)') 'Error while reading SC interaction',
1518 & 'potential file - unknown potential type.'
1520 call MPI_Finalize(Ierror)
1526 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1527 & ', exponents are ',expon,2*expon
1528 goto (10,20,30,30,40) ipot
1529 C----------------------- LJ potential ---------------------------------
1530 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1531 & (sigma0(i),i=1,ntyp)
1533 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1534 write (iout,'(a/)') 'The epsilon array:'
1535 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1536 write (iout,'(/a)') 'One-body parameters:'
1537 write (iout,'(a,4x,a)') 'residue','sigma'
1538 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1541 C----------------------- LJK potential --------------------------------
1542 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1543 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1545 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1546 write (iout,'(a/)') 'The epsilon array:'
1547 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1548 write (iout,'(/a)') 'One-body parameters:'
1549 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1550 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1554 C---------------------- GB or BP potential -----------------------------
1556 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1558 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1559 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1560 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1561 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1562 C now we start reading lipid
1564 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1565 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1567 C print *,"WARNING!!"
1569 C epslip(i,j)=epslip(i,j)+0.05d0
1572 if (lprint) write(iout,*) epslip(1,1),"OK?"
1573 C For the GB potential convert sigma'**2 into chi'
1576 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1580 write (iout,'(/a/)') 'Parameters of the BP potential:'
1581 write (iout,'(a/)') 'The epsilon array:'
1582 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1583 write (iout,'(/a)') 'One-body parameters:'
1584 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1586 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1587 & chip(i),alp(i),i=1,ntyp)
1590 C--------------------- GBV potential -----------------------------------
1591 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1592 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1593 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1595 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1596 write (iout,'(a/)') 'The epsilon array:'
1597 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1598 write (iout,'(/a)') 'One-body parameters:'
1599 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1600 & 's||/s_|_^2',' chip ',' alph '
1601 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1602 & sigii(i),chip(i),alp(i),i=1,ntyp)
1606 C-----------------------------------------------------------------------
1607 C Calculate the "working" parameters of SC interactions.
1611 epslip(i,j)=epslip(j,i)
1616 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1617 sigma(j,i)=sigma(i,j)
1618 rs0(i,j)=dwa16*sigma(i,j)
1622 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1623 & 'Working parameters of the SC interactions:',
1624 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1629 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1638 sigeps=dsign(1.0D0,epsij)
1640 aa_aq(i,j)=epsij*rrij*rrij
1641 bb_aq(i,j)=-sigeps*epsij*rrij
1642 aa_aq(j,i)=aa_aq(i,j)
1643 bb_aq(j,i)=bb_aq(i,j)
1644 epsijlip=epslip(i,j)
1645 sigeps=dsign(1.0D0,epsijlip)
1646 epsijlip=dabs(epsijlip)
1647 aa_lip(i,j)=epsijlip*rrij*rrij
1648 bb_lip(i,j)=-sigeps*epsijlip*rrij
1649 aa_lip(j,i)=aa_lip(i,j)
1650 bb_lip(j,i)=bb_lip(i,j)
1652 sigt1sq=sigma0(i)**2
1653 sigt2sq=sigma0(j)**2
1656 ratsig1=sigt2sq/sigt1sq
1657 ratsig2=1.0D0/ratsig1
1658 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1659 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1660 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1664 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1665 sigmaii(i,j)=rsum_max
1666 sigmaii(j,i)=rsum_max
1668 c sigmaii(i,j)=r0(i,j)
1669 c sigmaii(j,i)=r0(i,j)
1671 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1672 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1673 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1674 augm(i,j)=epsij*r_augm**(2*expon)
1675 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1682 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1683 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1684 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1689 write(iout,*) "tube param"
1690 read(itube,*) epspeptube,sigmapeptube
1691 sigmapeptube=sigmapeptube**6
1692 sigeps=dsign(1.0D0,epspeptube)
1693 epspeptube=dabs(epspeptube)
1694 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1695 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1696 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1698 read(itube,*) epssctube,sigmasctube
1699 sigmasctube=sigmasctube**6
1700 sigeps=dsign(1.0D0,epssctube)
1701 epssctube=dabs(epssctube)
1702 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1703 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1704 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1710 C Define the SC-p interaction constants (hard-coded; old style)
1713 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1715 c aad(i,1)=0.3D0*4.0D0**12
1716 C Following line for constants currently implemented
1717 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1718 aad(i,1)=1.5D0*4.0D0**12
1719 c aad(i,1)=0.17D0*5.6D0**12
1721 C "Soft" SC-p repulsion
1723 C Following line for constants currently implemented
1724 c aad(i,1)=0.3D0*4.0D0**6
1725 C "Hard" SC-p repulsion
1726 bad(i,1)=3.0D0*4.0D0**6
1727 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1736 C 8/9/01 Read the SC-p interaction constants from file
1739 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1742 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1743 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1744 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1745 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1749 write (iout,*) "Parameters of SC-p interactions:"
1751 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1752 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1758 C Define the constants of the disulfide bridge
1762 c Old arbitrary potential - commented out.
1767 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1768 c energy surface of diethyl disulfide.
1769 c A. Liwo and U. Kozlowska, 11/24/03
1785 C if(me.eq.king) then
1786 C write (iout,'(/a)') "Disulfide bridge parameters:"
1787 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1788 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1789 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1790 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1793 if (shield_mode.gt.0) then
1794 C VSolvSphere the volume of solving sphere
1795 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1796 C there will be no distinction between proline peptide group and normal peptide
1797 C group in case of shielding parameters
1798 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
1799 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1800 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1801 write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
1803 C long axis of side chain
1805 long_r_sidechain(i)=vbldsc0(1,i)
1806 short_r_sidechain(i)=sigma0(i)
1811 111 write (iout,*) "Error reading bending energy parameters."
1813 112 write (iout,*) "Error reading rotamer energy parameters."
1815 113 write (iout,*) "Error reading torsional energy parameters."
1817 114 write (iout,*) "Error reading double torsional energy parameters."
1820 & "Error reading cumulant (multibody energy) parameters."
1822 116 write (iout,*) "Error reading electrostatic energy parameters."
1824 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1826 117 write (iout,*) "Error reading side chain interaction parameters."
1828 118 write (iout,*) "Error reading SCp interaction parameters."
1830 119 write (iout,*) "Error reading SCCOR parameters"
1832 121 write (iout,*) "Error in Czybyshev parameters"
1835 call MPI_Finalize(Ierror)
1842 subroutine getenv_loc(var, val)
1843 character(*) var, val
1846 character(2000) line
1849 open (196,file='env',status='old',readonly,shared)
1851 c write(*,*)'looking for ',var
1852 10 read(196,*,err=11,end=11)line
1853 iread=index(line,var)
1854 c write(*,*)iread,' ',var,' ',line
1855 if (iread.eq.0) go to 10
1856 c write(*,*)'---> ',line
1862 iread=iread+ilen(var)+1
1863 read (line(iread:),*,err=12,end=12) val
1864 c write(*,*)'OK: ',var,' = ',val
1870 #elif (defined CRAY)
1871 integer lennam,lenval,ierror
1873 c getenv using a POSIX call, useful on the T3D
1874 c Sept 1996, comment out error check on advice of H. Pritchard
1877 if(lennam.le.0) stop '--error calling getenv--'
1878 call pxfgetenv(var,lennam,val,lenval,ierror)
1879 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1881 call getenv(var,val)
1883 C set the variables used for shielding effect
1884 C if (shield_mode.gt.0) then
1885 C VSolvSphere the volume of solving sphere
1887 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1888 C there will be no distinction between proline peptide group and normal peptide
1889 C group in case of shielding parameters
1890 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1891 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1892 C long axis of side chain
1894 C long_r_sidechain(i)=vbldsc0(1,i)
1895 C short_r_sidechain(i)=sigma0(i)
1897 C lets set the buffor value
1902 c---------------------------------------------------------------------------
1903 subroutine weightread
1905 C Read the parameters of the probability distributions of the virtual-bond
1906 C valence angles and the side chains and energy parameters.
1908 C Important! Energy-term weights ARE NOT read here; they are read from the
1909 C main input file instead, because NO defaults have yet been set for these
1913 include 'DIMENSIONS'
1918 include 'COMMON.IOUNITS'
1919 include 'COMMON.CHAIN'
1920 include 'COMMON.INTERACT'
1921 include 'COMMON.GEO'
1922 include 'COMMON.LOCAL'
1923 include 'COMMON.TORSION'
1924 include 'COMMON.SCCOR'
1925 include 'COMMON.SCROT'
1926 include 'COMMON.FFIELD'
1927 include 'COMMON.NAMES'
1928 include 'COMMON.SBRIDGE'
1930 include 'COMMON.SETUP'
1931 include 'COMMON.CONTROL'
1932 include 'COMMON.SHIELD'
1933 character*1000 weightcard
1935 double precision scalscp,wlong
1937 c READ energy-term weights
1939 call card_concat(weightcard)
1940 call reada(weightcard,'WLONG',wlong,1.0D0)
1941 call reada(weightcard,'WSC',wsc,wlong)
1942 call reada(weightcard,'WSCP',wscp,wlong)
1943 call reada(weightcard,'WELEC',welec,1.0D0)
1944 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1945 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1946 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1947 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1948 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1949 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1950 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1951 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1952 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1953 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1954 call reada(weightcard,'WBOND',wbond,1.0D0)
1955 call reada(weightcard,'WTOR',wtor,1.0D0)
1956 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1957 call reada(weightcard,'WANG',wang,1.0D0)
1958 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1959 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
1960 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
1961 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
1962 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
1963 call reada(weightcard,'SCAL14',scal14,0.4D0)
1964 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1965 call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
1966 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1967 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1968 call reada(weightcard,'TEMP0',temp0,300.0d0)
1969 call reada(weightcard,'WSHIELD',wshield,1.0d0)
1970 call reada(weightcard,'WLT',wliptran,0.0D0)
1971 call reada(weightcard,'WTUBE',wtube,1.0D0)
1972 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
1973 if (index(weightcard,'SOFT').gt.0) ipot=6
1974 C 12/1/95 Added weight for the multi-body term WCORR
1975 call reada(weightcard,'WCORRH',wcorr,1.0D0)
1976 if (wcorr4.gt.0.0d0) wcorr=wcorr4
1997 weights(28)=wdfa_dist
1998 weights(29)=wdfa_tor
1999 weights(30)=wdfa_nei
2000 weights(31)=wdfa_beta
2002 aad(i,1)=scalscp*aad(i,1)
2003 aad(i,2)=scalscp*aad(i,2)
2004 bad(i,1)=scalscp*bad(i,1)
2005 bad(i,2)=scalscp*bad(i,2)
2007 if(me.eq.king.or..not.out1file)
2008 & write (iout,100) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
2009 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
2011 100 format (/'Energy-term weights (unscaled):'//
2012 & 'WSCC= ',f10.6,' (SC-SC)'/
2013 & 'WSCP= ',f10.6,' (SC-p)'/
2014 & 'WELEC= ',f10.6,' (p-p electr)'/
2015 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
2016 & 'WBOND= ',f10.6,' (stretching)'/
2017 & 'WANG= ',f10.6,' (bending)'/
2018 & 'WSCLOC= ',f10.6,' (SC local)'/
2019 & 'WTOR= ',f10.6,' (torsional)'/
2020 & 'WTORD= ',f10.6,' (double torsional)'/
2021 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
2022 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
2023 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
2024 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
2025 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
2026 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
2027 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
2028 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
2029 & 'WTURN6= ',f10.6,' (turns, 6th order)')
2030 if(me.eq.king.or..not.out1file)then
2031 if (wcorr4.gt.0.0d0) then
2032 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
2033 & 'between contact pairs of peptide groups'
2034 write (iout,'(2(a,f5.3/))')
2035 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
2036 & 'Range of quenching the correlation terms:',2*delt_corr
2037 else if (wcorr.gt.0.0d0) then
2038 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
2039 & 'between contact pairs of peptide groups'
2041 write (iout,'(a,f8.3)')
2042 & 'Scaling factor of 1,4 SC-p interactions:',scal14
2043 write (iout,'(a,f8.3)')
2044 & 'General scaling factor of SC-p interactions:',scalscp
2046 r0_corr=cutoff_corr-delt_corr
2048 call rescale_weights(t_bath)
2049 if(me.eq.king.or..not.out1file)
2050 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
2051 & wtor_d,wstrain,wel_loc,
2053 & wcorr,wcorr5,wcorr6,
2060 22 format (/'Energy-term weights (scaled):'//
2061 & 'WSCC= ',f10.6,' (SC-SC)'/
2062 & 'WSCP= ',f10.6,' (SC-p)'/
2063 & 'WELEC= ',f10.6,' (p-p electr)'/
2064 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
2065 & 'WBOND= ',f10.6,' (stretching)'/
2066 & 'WANG= ',f10.6,' (bending)'/
2067 & 'WSCLOC= ',f10.6,' (SC local)'/
2068 & 'WTOR= ',f10.6,' (torsional)'/
2069 & 'WTORD= ',f10.6,' (double torsional)'/
2070 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
2071 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
2073 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
2074 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
2075 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
2077 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
2078 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
2079 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
2081 & 'WTURN6= ',f10.6,' (turns, 6th order)'
2084 if(me.eq.king.or..not.out1file)
2085 & write (iout,*) "Reference temperature for weights calculation:",
2088 write (iout,'(/a/)') "DFA pseudopotential parameters:"
2089 write (iout,'(a,f10.6,a)')
2090 & "WDFAD= ",wdfa_dist," (distance)",
2091 & "WDFAT= ",wdfa_tor," (backbone angles)",
2092 & "WDFAN= ",wdfa_nei," (neighbors)",
2093 & "WDFAB= ",wdfa_beta," (beta structure)"
2095 call reada(weightcard,"D0CM",d0cm,3.78d0)
2096 call reada(weightcard,"AKCM",akcm,15.1d0)
2097 call reada(weightcard,"AKTH",akth,11.0d0)
2098 call reada(weightcard,"AKCT",akct,12.0d0)
2099 call reada(weightcard,"V1SS",v1ss,-1.08d0)
2100 call reada(weightcard,"V2SS",v2ss,7.61d0)
2101 call reada(weightcard,"V3SS",v3ss,13.7d0)
2102 call reada(weightcard,"EBR",ebr,-5.50D0)
2103 call reada(weightcard,"ATRISS",atriss,0.301D0)
2104 call reada(weightcard,"BTRISS",btriss,0.021D0)
2105 call reada(weightcard,"CTRISS",ctriss,1.001D0)
2106 call reada(weightcard,"DTRISS",dtriss,1.001D0)
2107 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
2109 dyn_ss_mask(i)=.false.
2113 dyn_ssbond_ij(i,j)=1.0d300
2117 call reada(weightcard,"HT",Ht,0.0D0)
2119 ss_depth=ebr/wsc-0.25*eps(1,1)
2120 Ht=Ht/wsc-0.25*eps(1,1)
2121 akcm=akcm*wstrain/wsc
2122 akth=akth*wstrain/wsc
2123 akct=akct*wstrain/wsc
2124 v1ss=v1ss*wstrain/wsc
2125 v2ss=v2ss*wstrain/wsc
2126 v3ss=v3ss*wstrain/wsc
2128 if (wstrain.ne.0.0) then
2129 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
2134 write (iout,*) "wshield,", wshield
2135 if(me.eq.king.or..not.out1file) then
2136 write (iout,*) "Parameters of the SS-bond potential:"
2137 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
2139 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
2140 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
2141 write (iout,*)" HT",Ht
2142 write (iout,*) "Parameters of the 'trisulfide' potential"
2143 write (iout,*) "ATRISS=", atriss
2144 write (iout,*) "BTRISS=", btriss
2145 write (iout,*) "CTRISS=", ctriss
2146 write (iout,*) "DTRISS=", dtriss
2147 c print *,'indpdb=',indpdb,' pdbref=',pdbref