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
10 implicit real*8 (a-h,o-z)
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'
28 include 'COMMON.SETUP'
29 include 'COMMON.CONTROL'
30 include 'COMMON.SHIELD'
32 character*1 onelett(4) /"G","A","P","D"/
33 character*1 toronelet(-2:2) /"p","a","G","A","P"/
35 dimension blower(3,3,maxlob)
38 character*3 lancuch,ucase
39 character*1000 weightcard
41 C For printing parameters after they are read set the following in the UNRES
44 C setenv PRINT_PARM YES
46 C To print parameters in LaTeX format rather than as ASCII tables:
50 call getenv_loc("PRINT_PARM",lancuch)
51 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
52 call getenv_loc("LATEX",lancuch)
53 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
55 dwa16=2.0d0**(1.0d0/6.0d0)
58 C Assign virtual-bond length
63 c Read the virtual-bond parameters, masses, and moments of inertia
64 c and Stokes radii of the peptide group and side chains
67 read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
70 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
75 dsc_inv(i)=1.0D0/dsc(i)
79 read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
81 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
82 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
87 dsc_inv(i)=1.0D0/dsc(i)
92 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
93 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
95 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
97 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
98 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
100 write (iout,'(13x,3f10.5)')
101 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
105 C reading lipid parameters
107 write (iout,*) "iliptranpar",iliptranpar
110 read(iliptranpar,*) pepliptran
112 read(iliptranpar,*) liptranene(i)
117 C Read the parameters of the probability distribution/energy expression
118 C of the virtual-bond valence angles theta
121 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
122 & (bthet(j,i,1,1),j=1,2)
123 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
124 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
125 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
129 athet(1,i,1,-1)=athet(1,i,1,1)
130 athet(2,i,1,-1)=athet(2,i,1,1)
131 bthet(1,i,1,-1)=-bthet(1,i,1,1)
132 bthet(2,i,1,-1)=-bthet(2,i,1,1)
133 athet(1,i,-1,1)=-athet(1,i,1,1)
134 athet(2,i,-1,1)=-athet(2,i,1,1)
135 bthet(1,i,-1,1)=bthet(1,i,1,1)
136 bthet(2,i,-1,1)=bthet(2,i,1,1)
140 athet(1,i,-1,-1)=athet(1,-i,1,1)
141 athet(2,i,-1,-1)=-athet(2,-i,1,1)
142 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
143 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
144 athet(1,i,-1,1)=athet(1,-i,1,1)
145 athet(2,i,-1,1)=-athet(2,-i,1,1)
146 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
147 bthet(2,i,-1,1)=bthet(2,-i,1,1)
148 athet(1,i,1,-1)=-athet(1,-i,1,1)
149 athet(2,i,1,-1)=athet(2,-i,1,1)
150 bthet(1,i,1,-1)=bthet(1,-i,1,1)
151 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
156 polthet(j,i)=polthet(j,-i)
159 gthet(j,i)=gthet(j,-i)
167 & 'Parameters of the virtual-bond valence angles:'
168 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
169 & ' ATHETA0 ',' A1 ',' A2 ',
172 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
173 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
175 write (iout,'(/a/9x,5a/79(1h-))')
176 & 'Parameters of the expression for sigma(theta_c):',
177 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
178 & ' ALPH3 ',' SIGMA0C '
180 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
181 & (polthet(j,i),j=0,3),sigc0(i)
183 write (iout,'(/a/9x,5a/79(1h-))')
184 & 'Parameters of the second gaussian:',
185 & ' THETA0 ',' SIGMA0 ',' G1 ',
188 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
189 & sig0(i),(gthet(j,i),j=1,3)
193 & 'Parameters of the virtual-bond valence angles:'
194 write (iout,'(/a/9x,5a/79(1h-))')
195 & 'Coefficients of expansion',
196 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
197 & ' b1*10^1 ',' b2*10^1 '
199 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
200 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
201 & (10*bthet(j,i,1,1),j=1,2)
203 write (iout,'(/a/9x,5a/79(1h-))')
204 & 'Parameters of the expression for sigma(theta_c):',
205 & ' alpha0 ',' alph1 ',' alph2 ',
206 & ' alhp3 ',' sigma0c '
208 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
209 & (polthet(j,i),j=0,3),sigc0(i)
211 write (iout,'(/a/9x,5a/79(1h-))')
212 & 'Parameters of the second gaussian:',
213 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
216 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
217 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
222 IF (tor_mode.eq.0) THEN
224 C Read the parameters of Utheta determined from ab initio surfaces
225 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
227 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
228 & ntheterm3,nsingle,ndouble
229 write (iout,*) "ithep",ithep
231 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
232 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
234 ithetyp(i)=-ithetyp(-i)
237 do i=-maxthetyp,maxthetyp
238 do j=-maxthetyp,maxthetyp
239 do k=-maxthetyp,maxthetyp
240 aa0thet(i,j,k,iblock)=0.0d0
242 aathet(l,i,j,k,iblock)=0.0d0
246 bbthet(m,l,i,j,k,iblock)=0.0d0
247 ccthet(m,l,i,j,k,iblock)=0.0d0
248 ddthet(m,l,i,j,k,iblock)=0.0d0
249 eethet(m,l,i,j,k,iblock)=0.0d0
255 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
256 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
264 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
266 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
267 c VAR:1=non-glicyne non-proline 2=proline
268 c VAR:negative values for D-aminoacid
270 do j=-nthetyp,nthetyp
271 do k=-nthetyp,nthetyp
272 read (ithep,'(6a)',end=111,err=111) res1
273 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
274 c VAR: aa0thet is variable describing the average value of Foureir
275 c VAR: expansion series
276 c VAR: aathet is foureir expansion in theta/2 angle for full formula
277 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
278 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
279 read (ithep,*,end=111,err=111)
280 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
281 read (ithep,*,end=111,err=111)
282 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
283 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
284 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
285 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
287 read (ithep,*,end=111,err=111)
288 & (((ffthet(llll,lll,ll,i,j,k,iblock),
289 & ffthet(lll,llll,ll,i,j,k,iblock),
290 & ggthet(llll,lll,ll,i,j,k,iblock),
291 & ggthet(lll,llll,ll,i,j,k,iblock),
292 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
297 C For dummy ends assign glycine-type coefficients of theta-only terms; the
298 C coefficients of theta-and-gamma-dependent terms are zero.
299 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
300 C RECOMENTDED AFTER VERSION 3.3)
304 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
305 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
307 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
308 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
311 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
313 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
316 C AND COMMENT THE LOOPS BELOW
320 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
321 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
323 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
324 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
327 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
329 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
333 C Substitution for D aminoacids from symmetry.
336 do j=-nthetyp,nthetyp
337 do k=-nthetyp,nthetyp
338 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
340 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
344 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
345 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
346 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
347 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
353 ffthet(llll,lll,ll,i,j,k,iblock)=
354 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
355 ffthet(lll,llll,ll,i,j,k,iblock)=
356 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
357 ggthet(llll,lll,ll,i,j,k,iblock)=
358 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
359 ggthet(lll,llll,ll,i,j,k,iblock)=
360 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
369 C Control printout of the coefficients of virtual-bond-angle potentials
372 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
375 do j=-nthetyp,nthetyp
376 do k=-nthetyp,nthetyp
377 write (iout,'(//4a)')
378 & 'Type ',toronelet(i),toronelet(j),toronelet(k)
379 write (iout,'(//a,10x,a)') " l","a[l]"
380 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
381 write (iout,'(i2,1pe15.5)')
382 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
384 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
385 & "b",l,"c",l,"d",l,"e",l
387 write (iout,'(i2,4(1pe15.5))') m,
388 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
389 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
393 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
394 & "f+",l,"f-",l,"g+",l,"g-",l
397 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
398 & ffthet(n,m,l,i,j,k,iblock),
399 & ffthet(m,n,l,i,j,k,iblock),
400 & ggthet(n,m,l,i,j,k,iblock),
401 & ggthet(m,n,l,i,j,k,iblock)
414 C here will be the apropriate recalibrating for D-aminoacid
415 read (ithep,*,end=121,err=121) nthetyp
416 do i=-nthetyp+1,nthetyp-1
417 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
418 do j=0,nbend_kcc_Tb(i)
419 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
424 & "Parameters of the valence-only potentials"
425 do i=-nthetyp+1,nthetyp-1
426 write (iout,'(2a)') "Type ",toronelet(i)
427 do k=0,nbend_kcc_Tb(i)
428 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
435 c write (2,*) "Start reading THETA_PDB",ithep_pdb
438 read (ithep_pdb,*,err=111,end=111)
439 & a0thet(i),(athet(j,i,1,1),j=1,2),
440 & (bthet(j,i,1,1),j=1,2)
441 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
442 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
443 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
447 athet(1,i,1,-1)=athet(1,i,1,1)
448 athet(2,i,1,-1)=athet(2,i,1,1)
449 bthet(1,i,1,-1)=-bthet(1,i,1,1)
450 bthet(2,i,1,-1)=-bthet(2,i,1,1)
451 athet(1,i,-1,1)=-athet(1,i,1,1)
452 athet(2,i,-1,1)=-athet(2,i,1,1)
453 bthet(1,i,-1,1)=bthet(1,i,1,1)
454 bthet(2,i,-1,1)=bthet(2,i,1,1)
458 athet(1,i,-1,-1)=athet(1,-i,1,1)
459 athet(2,i,-1,-1)=-athet(2,-i,1,1)
460 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
461 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
462 athet(1,i,-1,1)=athet(1,-i,1,1)
463 athet(2,i,-1,1)=-athet(2,-i,1,1)
464 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
465 bthet(2,i,-1,1)=bthet(2,-i,1,1)
466 athet(1,i,1,-1)=-athet(1,-i,1,1)
467 athet(2,i,1,-1)=athet(2,-i,1,1)
468 bthet(1,i,1,-1)=bthet(1,-i,1,1)
469 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
474 polthet(j,i)=polthet(j,-i)
477 gthet(j,i)=gthet(j,-i)
480 c write (2,*) "End reading THETA_PDB"
486 C Read the parameters of the probability distribution/energy expression
487 C of the side chains.
490 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
494 dsc_inv(i)=1.0D0/dsc(i)
505 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
506 & ((blower(k,l,1),l=1,k),k=1,3)
507 censc(1,1,-i)=censc(1,1,i)
508 censc(2,1,-i)=censc(2,1,i)
509 censc(3,1,-i)=-censc(3,1,i)
511 read (irotam,*,end=112,err=112) bsc(j,i)
512 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
513 & ((blower(k,l,j),l=1,k),k=1,3)
514 censc(1,j,-i)=censc(1,j,i)
515 censc(2,j,-i)=censc(2,j,i)
516 censc(3,j,-i)=-censc(3,j,i)
517 C BSC is amplitude of Gaussian
524 akl=akl+blower(k,m,j)*blower(l,m,j)
528 if (((k.eq.3).and.(l.ne.3))
529 & .or.((l.eq.3).and.(k.ne.3))) then
530 gaussc(k,l,j,-i)=-akl
531 gaussc(l,k,j,-i)=-akl
543 write (iout,'(/a)') 'Parameters of side-chain local geometry'
548 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
549 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
550 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
551 & 'log h',(bsc(j,i),j=1,nlobi)
552 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
553 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
555 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
556 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
559 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
560 write (iout,'(a,f10.4,4(16x,f10.4))')
561 & 'Center ',(bsc(j,i),j=1,nlobi)
562 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
571 C Read scrot parameters for potentials determined from all-atom AM1 calculations
572 C added by Urszula Kozlowska 07/11/2007
575 read (irotam,*,end=112,err=112)
577 read (irotam,*,end=112,err=112)
580 read(irotam,*,end=112,err=112) sc_parmin(j,i)
585 C Read the parameters of the probability distribution/energy expression
586 C of the side chains.
588 write (2,*) "Start reading ROTAM_PDB"
590 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
594 dsc_inv(i)=1.0D0/dsc(i)
605 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
606 & ((blower(k,l,1),l=1,k),k=1,3)
608 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
609 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
610 & ((blower(k,l,j),l=1,k),k=1,3)
617 akl=akl+blower(k,m,j)*blower(l,m,j)
627 write (2,*) "End reading ROTAM_PDB"
631 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
632 C interaction energy of the Gly, Ala, and Pro prototypes.
634 read (ifourier,*) nloctyp
635 SPLIT_FOURIERTOR = nloctyp.lt.0
636 nloctyp = iabs(nloctyp)
638 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
639 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
640 itype2loc(ntyp1)=nloctyp
641 iloctyp(nloctyp)=ntyp1
643 itype2loc(-i)=-itype2loc(i)
652 iloctyp(-i)=-iloctyp(i)
654 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
655 c write (iout,*) "nloctyp",nloctyp,
656 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
659 c write (iout,*) "NEWCORR",i
660 read (ifourier,*,end=115,err=115)
663 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
666 c write (iout,*) "NEWCORR BNEW1"
667 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
670 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
673 c write (iout,*) "NEWCORR BNEW2"
674 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
676 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
677 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
679 c write (iout,*) "NEWCORR CCNEW"
680 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
682 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
683 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
685 c write (iout,*) "NEWCORR DDNEW"
686 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
690 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
694 c write (iout,*) "NEWCORR EENEW1"
695 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
697 read (ifourier,*,end=115,err=115) e0new(ii,i)
699 c write (iout,*) (e0new(ii,i),ii=1,3)
701 c write (iout,*) "NEWCORR EENEW"
704 ccnew(ii,1,i)=ccnew(ii,1,i)/2
705 ccnew(ii,2,i)=ccnew(ii,2,i)/2
706 ddnew(ii,1,i)=ddnew(ii,1,i)/2
707 ddnew(ii,2,i)=ddnew(ii,2,i)/2
712 bnew1(ii,1,-i)= bnew1(ii,1,i)
713 bnew1(ii,2,-i)=-bnew1(ii,2,i)
714 bnew2(ii,1,-i)= bnew2(ii,1,i)
715 bnew2(ii,2,-i)=-bnew2(ii,2,i)
718 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
719 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
720 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
721 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
722 ccnew(ii,1,-i)=ccnew(ii,1,i)
723 ccnew(ii,2,-i)=-ccnew(ii,2,i)
724 ddnew(ii,1,-i)=ddnew(ii,1,i)
725 ddnew(ii,2,-i)=-ddnew(ii,2,i)
727 e0new(1,-i)= e0new(1,i)
728 e0new(2,-i)=-e0new(2,i)
729 e0new(3,-i)=-e0new(3,i)
731 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
732 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
733 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
734 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
738 write (iout,'(a)') "Coefficients of the multibody terms"
739 do i=-nloctyp+1,nloctyp-1
740 write (iout,*) "Type: ",onelet(iloctyp(i))
741 write (iout,*) "Coefficients of the expansion of B1"
743 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
745 write (iout,*) "Coefficients of the expansion of B2"
747 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
749 write (iout,*) "Coefficients of the expansion of C"
750 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
751 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
752 write (iout,*) "Coefficients of the expansion of D"
753 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
754 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
755 write (iout,*) "Coefficients of the expansion of E"
756 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
759 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
764 IF (SPLIT_FOURIERTOR) THEN
766 c write (iout,*) "NEWCORR TOR",i
767 read (ifourier,*,end=115,err=115)
770 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
773 c write (iout,*) "NEWCORR BNEW1 TOR"
774 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
777 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
780 c write (iout,*) "NEWCORR BNEW2 TOR"
781 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
783 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
784 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
786 c write (iout,*) "NEWCORR CCNEW TOR"
787 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
789 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
790 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
792 c write (iout,*) "NEWCORR DDNEW TOR"
793 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
797 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
801 c write (iout,*) "NEWCORR EENEW1 TOR"
802 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
804 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
806 c write (iout,*) (e0newtor(ii,i),ii=1,3)
808 c write (iout,*) "NEWCORR EENEW TOR"
811 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
812 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
813 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
814 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
819 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
820 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
821 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
822 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
825 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
826 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
827 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
828 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
830 e0newtor(1,-i)= e0newtor(1,i)
831 e0newtor(2,-i)=-e0newtor(2,i)
832 e0newtor(3,-i)=-e0newtor(3,i)
834 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
835 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
836 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
837 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
842 & "Single-body coefficients of the torsional potentials"
843 do i=-nloctyp+1,nloctyp-1
844 write (iout,*) "Type: ",onelet(iloctyp(i))
845 write (iout,*) "Coefficients of the expansion of B1tor"
847 write (iout,'(3hB1(,i1,1h),3f10.5)')
848 & j,(bnew1tor(k,j,i),k=1,3)
850 write (iout,*) "Coefficients of the expansion of B2tor"
852 write (iout,'(3hB2(,i1,1h),3f10.5)')
853 & j,(bnew2tor(k,j,i),k=1,3)
855 write (iout,*) "Coefficients of the expansion of Ctor"
856 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
857 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
858 write (iout,*) "Coefficients of the expansion of Dtor"
859 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
860 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
861 write (iout,*) "Coefficients of the expansion of Etor"
862 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
865 write (iout,'(1hE,2i1,2f10.5)')
866 & j,k,(eenewtor(l,j,k,i),l=1,2)
872 do i=-nloctyp+1,nloctyp-1
875 bnew1tor(ii,j,i)=bnew1(ii,j,i)
880 bnew2tor(ii,j,i)=bnew2(ii,j,i)
884 ccnewtor(ii,1,i)=ccnew(ii,1,i)
885 ccnewtor(ii,2,i)=ccnew(ii,2,i)
886 ddnewtor(ii,1,i)=ddnew(ii,1,i)
887 ddnewtor(ii,2,i)=ddnew(ii,2,i)
890 ENDIF !SPLIT_FOURIER_TOR
893 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
895 read (ifourier,*,end=115,err=115)
896 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
898 write (iout,*) 'Type ',onelet(iloctyp(i))
899 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
913 c B1tilde(1,i) = b(3)
914 c B1tilde(2,i) =-b(5)
915 c B1tilde(1,-i) =-b(3)
916 c B1tilde(2,-i) =b(5)
923 cc B1tilde(1,i) = b(3,i)
924 cc B1tilde(2,i) =-b(5,i)
925 C B1tilde(1,-i) =-b(3,i)
926 C B1tilde(2,-i) =b(5,i)
927 cc b1tilde(1,i)=0.0d0
928 cc b1tilde(2,i)=0.0d0
940 CCold(1,1,-i)= b(7,i)
941 CCold(2,2,-i)=-b(7,i)
942 CCold(2,1,-i)=-b(9,i)
943 CCold(1,2,-i)=-b(9,i)
948 c Ctilde(1,1,i)= CCold(1,1,i)
949 c Ctilde(1,2,i)= CCold(1,2,i)
950 c Ctilde(2,1,i)=-CCold(2,1,i)
951 c Ctilde(2,2,i)=-CCold(2,2,i)
953 c Ctilde(1,1,i)=0.0d0
954 c Ctilde(1,2,i)=0.0d0
955 c Ctilde(2,1,i)=0.0d0
956 c Ctilde(2,2,i)=0.0d0
961 DDold(1,1,-i)= b(6,i)
962 DDold(2,2,-i)=-b(6,i)
963 DDold(2,1,-i)=-b(8,i)
964 DDold(1,2,-i)=-b(8,i)
969 c Dtilde(1,1,i)= DD(1,1,i)
970 c Dtilde(1,2,i)= DD(1,2,i)
971 c Dtilde(2,1,i)=-DD(2,1,i)
972 c Dtilde(2,2,i)=-DD(2,2,i)
974 c Dtilde(1,1,i)=0.0d0
975 c Dtilde(1,2,i)=0.0d0
976 c Dtilde(2,1,i)=0.0d0
977 c Dtilde(2,2,i)=0.0d0
978 EEold(1,1,i)= b(10,i)+b(11,i)
979 EEold(2,2,i)=-b(10,i)+b(11,i)
980 EEold(2,1,i)= b(12,i)-b(13,i)
981 EEold(1,2,i)= b(12,i)+b(13,i)
982 EEold(1,1,-i)= b(10,i)+b(11,i)
983 EEold(2,2,-i)=-b(10,i)+b(11,i)
984 EEold(2,1,-i)=-b(12,i)+b(13,i)
985 EEold(1,2,-i)=-b(12,i)-b(13,i)
986 write(iout,*) "TU DOCHODZE"
992 c ee(2,1,i)=ee(1,2,i)
997 &"Coefficients of the cumulants (independent of valence angles)"
998 do i=-nloctyp+1,nloctyp-1
999 write (iout,*) 'Type ',onelet(iloctyp(i))
1001 write(iout,'(2f10.5)') B(3,i),B(5,i)
1003 write(iout,'(2f10.5)') B(2,i),B(4,i)
1006 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1010 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1014 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1022 C Read torsional parameters in old format
1024 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1025 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1026 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1031 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1037 write (iout,'(/a/)') 'Torsional constants:'
1040 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1041 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1047 C Read torsional parameters
1049 IF (TOR_MODE.eq.0) THEN
1051 read (itorp,*,end=113,err=113) ntortyp
1052 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1054 itype2loc(i)=itortyp(i)
1057 itype2loc(-i)=-itype2loc(i)
1059 itortyp(ntyp1)=ntortyp
1062 itortyp(i)=-itortyp(-i)
1064 write (iout,*) 'ntortyp',ntortyp
1066 do j=-ntortyp+1,ntortyp-1
1067 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1069 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1070 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1073 do k=1,nterm(i,j,iblock)
1074 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1076 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1077 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1078 v0ij=v0ij+si*v1(k,i,j,iblock)
1080 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
1081 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
1082 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
1084 do k=1,nlor(i,j,iblock)
1085 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1086 & vlor2(k,i,j),vlor3(k,i,j)
1087 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1090 v0(-i,-j,iblock)=v0ij
1096 write (iout,'(/a/)') 'Parameters of the SCCOR potentials:'
1099 do j=-ntortyp+1,ntortyp-1
1100 write (iout,*) 'ityp',i,' jtyp',j
1101 write (iout,*) 'Fourier constants'
1102 do k=1,nterm(i,j,iblock)
1103 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1106 write (iout,*) 'Lorenz constants'
1107 do k=1,nlor(i,j,iblock)
1108 write (iout,'(3(1pe15.5))')
1109 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1117 C 6/23/01 Read parameters for double torsionals
1121 do j=-ntortyp+1,ntortyp-1
1122 do k=-ntortyp+1,ntortyp-1
1123 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1124 c write (iout,*) "OK onelett",
1127 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1128 & .or. t3.ne.toronelet(k)) then
1129 write (iout,*) "Error in double torsional parameter file",
1132 call MPI_Finalize(Ierror)
1134 stop "Error in double torsional parameter file"
1136 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1137 & ntermd_2(i,j,k,iblock)
1138 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1139 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1140 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1141 & ntermd_1(i,j,k,iblock))
1142 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1143 & ntermd_1(i,j,k,iblock))
1144 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1145 & ntermd_1(i,j,k,iblock))
1146 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1147 & ntermd_1(i,j,k,iblock))
1148 C Martix of D parameters for one dimesional foureir series
1149 do l=1,ntermd_1(i,j,k,iblock)
1150 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1151 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1152 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1153 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1154 c write(iout,*) "whcodze" ,
1155 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1157 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1158 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1159 & v2s(m,l,i,j,k,iblock),
1160 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1161 C Martix of D parameters for two dimesional fourier series
1162 do l=1,ntermd_2(i,j,k,iblock)
1164 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1165 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1166 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1167 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1176 write (iout,*) 'Constants for double torsionals'
1179 do j=-ntortyp+1,ntortyp-1
1180 do k=-ntortyp+1,ntortyp-1
1181 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1182 & ' nsingle',ntermd_1(i,j,k,iblock),
1183 & ' ndouble',ntermd_2(i,j,k,iblock)
1185 write (iout,*) 'Single angles:'
1186 do l=1,ntermd_1(i,j,k,iblock)
1187 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1188 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1189 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1190 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1193 write (iout,*) 'Pairs of angles:'
1194 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1195 do l=1,ntermd_2(i,j,k,iblock)
1196 write (iout,'(i5,20f10.5)')
1197 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1200 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1201 do l=1,ntermd_2(i,j,k,iblock)
1202 write (iout,'(i5,20f10.5)')
1203 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1204 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1213 ELSE IF (TOR_MODE.eq.1) THEN
1215 C read valence-torsional parameters
1216 read (itorp,*,end=121,err=121) ntortyp
1218 write (iout,*) "Valence-torsional parameters read in ntortyp",
1220 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
1221 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1224 itype2loc(i)=itortyp(i)
1228 itortyp(i)=-itortyp(-i)
1230 do i=-ntortyp+1,ntortyp-1
1231 do j=-ntortyp+1,ntortyp-1
1232 C first we read the cos and sin gamma parameters
1233 read (itorp,'(13x,a)',end=121,err=121) string
1234 write (iout,*) i,j,string
1235 read (itorp,*,end=121,err=121)
1236 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1237 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1238 do k=1,nterm_kcc(j,i)
1239 do l=1,nterm_kcc_Tb(j,i)
1240 do ll=1,nterm_kcc_Tb(j,i)
1241 read (itorp,*,end=121,err=121) ii,jj,kk,
1242 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1250 c AL 4/8/16: Calculate coefficients from one-body parameters
1253 itortyp(i)=itype2loc(i)
1256 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1258 do i=-ntortyp+1,ntortyp-1
1259 do j=-ntortyp+1,ntortyp-1
1262 do k=1,nterm_kcc_Tb(j,i)
1263 do l=1,nterm_kcc_Tb(j,i)
1264 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1265 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1266 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1267 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1270 do k=1,nterm_kcc_Tb(j,i)
1271 do l=1,nterm_kcc_Tb(j,i)
1273 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1274 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1275 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1276 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1278 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1279 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1280 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1281 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1285 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)
1289 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1294 if (tor_mode.gt.0 .and. lprint) then
1295 c Print valence-torsional parameters
1297 & "Parameters of the valence-torsional potentials"
1298 do i=-ntortyp+1,ntortyp-1
1299 do j=-ntortyp+1,ntortyp-1
1300 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1301 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1302 do k=1,nterm_kcc(j,i)
1303 do l=1,nterm_kcc_Tb(j,i)
1304 do ll=1,nterm_kcc_Tb(j,i)
1305 write (iout,'(3i5,2f15.4)')
1306 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1315 C Read of Side-chain backbone correlation parameters
1316 C Modified 11 May 2012 by Adasko
1319 read (isccor,*,end=119,err=119) nsccortyp
1321 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1323 isccortyp(i)=-isccortyp(-i)
1325 iscprol=isccortyp(20)
1326 c write (iout,*) 'ntortyp',ntortyp
1328 cc maxinter is maximum interaction sites
1332 read (isccor,*,end=119,err=119)
1333 &nterm_sccor(i,j),nlor_sccor(i,j)
1339 nterm_sccor(-i,j)=nterm_sccor(i,j)
1340 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1341 nterm_sccor(i,-j)=nterm_sccor(i,j)
1342 do k=1,nterm_sccor(i,j)
1343 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1345 if (j.eq.iscprol) then
1346 if (i.eq.isccortyp(10)) then
1347 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1348 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1350 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1351 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1352 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1353 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1354 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1355 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1356 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1357 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1360 if (i.eq.isccortyp(10)) then
1361 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1362 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1364 if (j.eq.isccortyp(10)) then
1365 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1366 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1368 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1369 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1370 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1371 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1372 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1373 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1377 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1378 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1379 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1380 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1383 do k=1,nlor_sccor(i,j)
1384 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1385 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1386 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1387 &(1+vlor3sccor(k,i,j)**2)
1389 v0sccor(l,i,j)=v0ijsccor
1390 v0sccor(l,-i,j)=v0ijsccor1
1391 v0sccor(l,i,-j)=v0ijsccor2
1392 v0sccor(l,-i,-j)=v0ijsccor3
1398 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1399 c write (iout,*) 'ntortyp',ntortyp
1401 cc maxinter is maximum interaction sites
1405 read (isccor,*,end=119,err=119)
1406 & nterm_sccor(i,j),nlor_sccor(i,j)
1410 do k=1,nterm_sccor(i,j)
1411 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1413 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1416 do k=1,nlor_sccor(i,j)
1417 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1418 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1419 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1420 &(1+vlor3sccor(k,i,j)**2)
1422 v0sccor(l,i,j)=v0ijsccor
1430 write (iout,'(/a/)') 'Torsional constants:'
1434 write (iout,*) 'ityp',i,' jtyp',j
1435 write (iout,*) 'Fourier constants'
1436 do k=1,nterm_sccor(i,j)
1437 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1439 write (iout,*) 'Lorenz constants'
1440 do k=1,nlor_sccor(i,j)
1441 write (iout,'(3(1pe15.5))')
1442 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1450 C Read electrostatic-interaction parameters
1454 write (iout,'(/a)') 'Electrostatic interaction constants:'
1455 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1456 & 'IT','JT','APP','BPP','AEL6','AEL3'
1458 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1459 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1460 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1461 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1466 app (i,j)=epp(i,j)*rri*rri
1467 bpp (i,j)=-2.0D0*epp(i,j)*rri
1468 ael6(i,j)=elpp6(i,j)*4.2D0**6
1469 ael3(i,j)=elpp3(i,j)*4.2D0**3
1471 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1472 & ael6(i,j),ael3(i,j)
1477 C Read side-chain interaction parameters.
1479 read (isidep,*,end=117,err=117) ipot,expon
1480 if (ipot.lt.1 .or. ipot.gt.5) then
1481 write (iout,'(2a)') 'Error while reading SC interaction',
1482 & 'potential file - unknown potential type.'
1484 call MPI_Finalize(Ierror)
1490 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1491 & ', exponents are ',expon,2*expon
1492 goto (10,20,30,30,40) ipot
1493 C----------------------- LJ potential ---------------------------------
1494 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1495 & (sigma0(i),i=1,ntyp)
1497 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1498 write (iout,'(a/)') 'The epsilon array:'
1499 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1500 write (iout,'(/a)') 'One-body parameters:'
1501 write (iout,'(a,4x,a)') 'residue','sigma'
1502 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1505 C----------------------- LJK potential --------------------------------
1506 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1507 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1509 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1510 write (iout,'(a/)') 'The epsilon array:'
1511 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1512 write (iout,'(/a)') 'One-body parameters:'
1513 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1514 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1518 C---------------------- GB or BP potential -----------------------------
1520 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1522 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1523 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1524 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1525 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1526 C now we start reading lipid
1528 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1529 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1531 C print *,"WARNING!!"
1533 C epslip(i,j)=epslip(i,j)+0.05d0
1536 if (lprint) write(iout,*) epslip(1,1),"OK?"
1537 C For the GB potential convert sigma'**2 into chi'
1540 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1544 write (iout,'(/a/)') 'Parameters of the BP potential:'
1545 write (iout,'(a/)') 'The epsilon array:'
1546 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1547 write (iout,'(/a)') 'One-body parameters:'
1548 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1550 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1551 & chip(i),alp(i),i=1,ntyp)
1554 C--------------------- GBV potential -----------------------------------
1555 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1556 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1557 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1559 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1560 write (iout,'(a/)') 'The epsilon array:'
1561 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1562 write (iout,'(/a)') 'One-body parameters:'
1563 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1564 & 's||/s_|_^2',' chip ',' alph '
1565 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1566 & sigii(i),chip(i),alp(i),i=1,ntyp)
1570 C-----------------------------------------------------------------------
1571 C Calculate the "working" parameters of SC interactions.
1575 epslip(i,j)=epslip(j,i)
1580 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1581 sigma(j,i)=sigma(i,j)
1582 rs0(i,j)=dwa16*sigma(i,j)
1586 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1587 & 'Working parameters of the SC interactions:',
1588 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1593 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1602 sigeps=dsign(1.0D0,epsij)
1604 aa_aq(i,j)=epsij*rrij*rrij
1605 bb_aq(i,j)=-sigeps*epsij*rrij
1606 aa_aq(j,i)=aa_aq(i,j)
1607 bb_aq(j,i)=bb_aq(i,j)
1608 epsijlip=epslip(i,j)
1609 sigeps=dsign(1.0D0,epsijlip)
1610 epsijlip=dabs(epsijlip)
1611 aa_lip(i,j)=epsijlip*rrij*rrij
1612 bb_lip(i,j)=-sigeps*epsijlip*rrij
1613 aa_lip(j,i)=aa_lip(i,j)
1614 bb_lip(j,i)=bb_lip(i,j)
1616 sigt1sq=sigma0(i)**2
1617 sigt2sq=sigma0(j)**2
1620 ratsig1=sigt2sq/sigt1sq
1621 ratsig2=1.0D0/ratsig1
1622 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1623 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1624 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1628 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1629 sigmaii(i,j)=rsum_max
1630 sigmaii(j,i)=rsum_max
1632 c sigmaii(i,j)=r0(i,j)
1633 c sigmaii(j,i)=r0(i,j)
1635 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1636 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1637 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1638 augm(i,j)=epsij*r_augm**(2*expon)
1639 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1646 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1647 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1648 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1653 write(iout,*) "tube param"
1654 read(itube,*) epspeptube,sigmapeptube
1655 sigmapeptube=sigmapeptube**6
1656 sigeps=dsign(1.0D0,epspeptube)
1657 epspeptube=dabs(epspeptube)
1658 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1659 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1660 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1662 read(itube,*) epssctube,sigmasctube
1663 sigmasctube=sigmasctube**6
1664 sigeps=dsign(1.0D0,epssctube)
1665 epssctube=dabs(epssctube)
1666 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1667 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1668 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1674 C Define the SC-p interaction constants (hard-coded; old style)
1677 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1679 c aad(i,1)=0.3D0*4.0D0**12
1680 C Following line for constants currently implemented
1681 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1682 aad(i,1)=1.5D0*4.0D0**12
1683 c aad(i,1)=0.17D0*5.6D0**12
1685 C "Soft" SC-p repulsion
1687 C Following line for constants currently implemented
1688 c aad(i,1)=0.3D0*4.0D0**6
1689 C "Hard" SC-p repulsion
1690 bad(i,1)=3.0D0*4.0D0**6
1691 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1700 C 8/9/01 Read the SC-p interaction constants from file
1703 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1706 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1707 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1708 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1709 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1713 write (iout,*) "Parameters of SC-p interactions:"
1715 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1716 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1722 C Define the constants of the disulfide bridge
1726 c Old arbitrary potential - commented out.
1731 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1732 c energy surface of diethyl disulfide.
1733 c A. Liwo and U. Kozlowska, 11/24/03
1749 C if(me.eq.king) then
1750 C write (iout,'(/a)') "Disulfide bridge parameters:"
1751 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1752 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1753 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1754 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1757 if (shield_mode.gt.0) then
1758 C VSolvSphere the volume of solving sphere
1759 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1760 C there will be no distinction between proline peptide group and normal peptide
1761 C group in case of shielding parameters
1762 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
1763 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1764 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1765 write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
1767 C long axis of side chain
1769 long_r_sidechain(i)=vbldsc0(1,i)
1770 short_r_sidechain(i)=sigma0(i)
1775 111 write (iout,*) "Error reading bending energy parameters."
1777 112 write (iout,*) "Error reading rotamer energy parameters."
1779 113 write (iout,*) "Error reading torsional energy parameters."
1781 114 write (iout,*) "Error reading double torsional energy parameters."
1784 & "Error reading cumulant (multibody energy) parameters."
1786 116 write (iout,*) "Error reading electrostatic energy parameters."
1788 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1790 117 write (iout,*) "Error reading side chain interaction parameters."
1792 118 write (iout,*) "Error reading SCp interaction parameters."
1794 119 write (iout,*) "Error reading SCCOR parameters"
1796 121 write (iout,*) "Error in Czybyshev parameters"
1799 call MPI_Finalize(Ierror)
1806 subroutine getenv_loc(var, val)
1807 character(*) var, val
1810 character(2000) line
1813 open (196,file='env',status='old',readonly,shared)
1815 c write(*,*)'looking for ',var
1816 10 read(196,*,err=11,end=11)line
1817 iread=index(line,var)
1818 c write(*,*)iread,' ',var,' ',line
1819 if (iread.eq.0) go to 10
1820 c write(*,*)'---> ',line
1826 iread=iread+ilen(var)+1
1827 read (line(iread:),*,err=12,end=12) val
1828 c write(*,*)'OK: ',var,' = ',val
1834 #elif (defined CRAY)
1835 integer lennam,lenval,ierror
1837 c getenv using a POSIX call, useful on the T3D
1838 c Sept 1996, comment out error check on advice of H. Pritchard
1841 if(lennam.le.0) stop '--error calling getenv--'
1842 call pxfgetenv(var,lennam,val,lenval,ierror)
1843 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1845 call getenv(var,val)
1847 C set the variables used for shielding effect
1848 C if (shield_mode.gt.0) then
1849 C VSolvSphere the volume of solving sphere
1851 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1852 C there will be no distinction between proline peptide group and normal peptide
1853 C group in case of shielding parameters
1854 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1855 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1856 C long axis of side chain
1858 C long_r_sidechain(i)=vbldsc0(1,i)
1859 C short_r_sidechain(i)=sigma0(i)
1861 C lets set the buffor value
1866 c---------------------------------------------------------------------------
1867 subroutine weightread
1869 C Read the parameters of the probability distributions of the virtual-bond
1870 C valence angles and the side chains and energy parameters.
1872 C Important! Energy-term weights ARE NOT read here; they are read from the
1873 C main input file instead, because NO defaults have yet been set for these
1876 implicit real*8 (a-h,o-z)
1877 include 'DIMENSIONS'
1882 include 'COMMON.IOUNITS'
1883 include 'COMMON.CHAIN'
1884 include 'COMMON.INTERACT'
1885 include 'COMMON.GEO'
1886 include 'COMMON.LOCAL'
1887 include 'COMMON.TORSION'
1888 include 'COMMON.SCCOR'
1889 include 'COMMON.SCROT'
1890 include 'COMMON.FFIELD'
1891 include 'COMMON.NAMES'
1892 include 'COMMON.SBRIDGE'
1894 include 'COMMON.SETUP'
1895 include 'COMMON.CONTROL'
1896 include 'COMMON.SHIELD'
1897 character*1000 weightcard
1899 c READ energy-term weights
1901 call card_concat(weightcard)
1902 call reada(weightcard,'WLONG',wlong,1.0D0)
1903 call reada(weightcard,'WSC',wsc,wlong)
1904 call reada(weightcard,'WSCP',wscp,wlong)
1905 call reada(weightcard,'WELEC',welec,1.0D0)
1906 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1907 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1908 write (iout,*) "wel_loc",wel_loc
1909 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1910 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1911 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1912 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1913 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1914 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1915 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1916 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1917 call reada(weightcard,'WBOND',wbond,1.0D0)
1918 call reada(weightcard,'WTOR',wtor,1.0D0)
1919 write (iout,*) "wtor",wtor
1920 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1921 call reada(weightcard,'WANG',wang,1.0D0)
1922 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1923 call reada(weightcard,'SCAL14',scal14,0.4D0)
1924 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1925 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1926 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1927 call reada(weightcard,'TEMP0',temp0,300.0d0)
1928 call reada(weightcard,'WSHIELD',wshield,1.0d0)
1929 call reada(weightcard,'WLT',wliptran,0.0D0)
1930 call reada(weightcard,'WTUBE',wtube,1.0D0)
1931 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
1932 if (index(weightcard,'SOFT').gt.0) ipot=6
1933 C 12/1/95 Added weight for the multi-body term WCORR
1934 call reada(weightcard,'WCORRH',wcorr,1.0D0)
1935 if (wcorr4.gt.0.0d0) wcorr=wcorr4
1957 aad(i,1)=scalscp*aad(i,1)
1958 aad(i,2)=scalscp*aad(i,2)
1959 bad(i,1)=scalscp*bad(i,1)
1960 bad(i,2)=scalscp*bad(i,2)
1962 if(me.eq.king.or..not.out1file)
1963 & write (iout,100) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
1964 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
1966 100 format (/'Energy-term weights (unscaled):'//
1967 & 'WSCC= ',f10.6,' (SC-SC)'/
1968 & 'WSCP= ',f10.6,' (SC-p)'/
1969 & 'WELEC= ',f10.6,' (p-p electr)'/
1970 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
1971 & 'WBOND= ',f10.6,' (stretching)'/
1972 & 'WANG= ',f10.6,' (bending)'/
1973 & 'WSCLOC= ',f10.6,' (SC local)'/
1974 & 'WTOR= ',f10.6,' (torsional)'/
1975 & 'WTORD= ',f10.6,' (double torsional)'/
1976 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
1977 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
1978 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
1979 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
1980 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
1981 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
1982 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
1983 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
1984 & 'WTURN6= ',f10.6,' (turns, 6th order)')
1985 if(me.eq.king.or..not.out1file)then
1986 if (wcorr4.gt.0.0d0) then
1987 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
1988 & 'between contact pairs of peptide groups'
1989 write (iout,'(2(a,f5.3/))')
1990 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
1991 & 'Range of quenching the correlation terms:',2*delt_corr
1992 else if (wcorr.gt.0.0d0) then
1993 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
1994 & 'between contact pairs of peptide groups'
1996 write (iout,'(a,f8.3)')
1997 & 'Scaling factor of 1,4 SC-p interactions:',scal14
1998 write (iout,'(a,f8.3)')
1999 & 'General scaling factor of SC-p interactions:',scalscp
2001 r0_corr=cutoff_corr-delt_corr
2003 call rescale_weights(t_bath)
2004 if(me.eq.king.or..not.out1file)
2005 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
2006 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
2008 22 format (/'Energy-term weights (scaled):'//
2009 & 'WSCC= ',f10.6,' (SC-SC)'/
2010 & 'WSCP= ',f10.6,' (SC-p)'/
2011 & 'WELEC= ',f10.6,' (p-p electr)'/
2012 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
2013 & 'WBOND= ',f10.6,' (stretching)'/
2014 & 'WANG= ',f10.6,' (bending)'/
2015 & 'WSCLOC= ',f10.6,' (SC local)'/
2016 & 'WTOR= ',f10.6,' (torsional)'/
2017 & 'WTORD= ',f10.6,' (double torsional)'/
2018 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
2019 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
2020 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
2021 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
2022 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
2023 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
2024 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
2025 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
2026 & 'WTURN6= ',f10.6,' (turns, 6th order)')
2027 if(me.eq.king.or..not.out1file)
2028 & write (iout,*) "Reference temperature for weights calculation:",
2030 call reada(weightcard,"D0CM",d0cm,3.78d0)
2031 call reada(weightcard,"AKCM",akcm,15.1d0)
2032 call reada(weightcard,"AKTH",akth,11.0d0)
2033 call reada(weightcard,"AKCT",akct,12.0d0)
2034 call reada(weightcard,"V1SS",v1ss,-1.08d0)
2035 call reada(weightcard,"V2SS",v2ss,7.61d0)
2036 call reada(weightcard,"V3SS",v3ss,13.7d0)
2037 call reada(weightcard,"EBR",ebr,-5.50D0)
2038 call reada(weightcard,"ATRISS",atriss,0.301D0)
2039 call reada(weightcard,"BTRISS",btriss,0.021D0)
2040 call reada(weightcard,"CTRISS",ctriss,1.001D0)
2041 call reada(weightcard,"DTRISS",dtriss,1.001D0)
2042 write (iout,*) "ATRISS=", atriss
2043 write (iout,*) "BTRISS=", btriss
2044 write (iout,*) "CTRISS=", ctriss
2045 write (iout,*) "DTRISS=", dtriss
2046 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
2048 dyn_ss_mask(i)=.false.
2052 dyn_ssbond_ij(i,j)=1.0d300
2055 call reada(weightcard,"HT",Ht,0.0D0)
2057 ss_depth=ebr/wsc-0.25*eps(1,1)
2058 Ht=Ht/wsc-0.25*eps(1,1)
2059 akcm=akcm*wstrain/wsc
2060 akth=akth*wstrain/wsc
2061 akct=akct*wstrain/wsc
2062 v1ss=v1ss*wstrain/wsc
2063 v2ss=v2ss*wstrain/wsc
2064 v3ss=v3ss*wstrain/wsc
2066 if (wstrain.ne.0.0) then
2067 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
2072 write (iout,*) "wshield,", wshield
2073 if(me.eq.king.or..not.out1file) then
2074 write (iout,*) "Parameters of the SS-bond potential:"
2075 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
2077 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
2078 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
2079 write (iout,*)" HT",Ht
2080 print *,'indpdb=',indpdb,' pdbref=',pdbref