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 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
230 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
232 ithetyp(i)=-ithetyp(-i)
235 do i=-maxthetyp,maxthetyp
236 do j=-maxthetyp,maxthetyp
237 do k=-maxthetyp,maxthetyp
238 aa0thet(i,j,k,iblock)=0.0d0
240 aathet(l,i,j,k,iblock)=0.0d0
244 bbthet(m,l,i,j,k,iblock)=0.0d0
245 ccthet(m,l,i,j,k,iblock)=0.0d0
246 ddthet(m,l,i,j,k,iblock)=0.0d0
247 eethet(m,l,i,j,k,iblock)=0.0d0
253 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
254 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
262 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
264 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
265 c VAR:1=non-glicyne non-proline 2=proline
266 c VAR:negative values for D-aminoacid
268 do j=-nthetyp,nthetyp
269 do k=-nthetyp,nthetyp
270 read (ithep,'(6a)',end=111,err=111) res1
271 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
272 c VAR: aa0thet is variable describing the average value of Foureir
273 c VAR: expansion series
274 c VAR: aathet is foureir expansion in theta/2 angle for full formula
275 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
276 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
277 read (ithep,*,end=111,err=111)
278 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
279 read (ithep,*,end=111,err=111)
280 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
281 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
282 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
283 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
285 read (ithep,*,end=111,err=111)
286 & (((ffthet(llll,lll,ll,i,j,k,iblock),
287 & ffthet(lll,llll,ll,i,j,k,iblock),
288 & ggthet(llll,lll,ll,i,j,k,iblock),
289 & ggthet(lll,llll,ll,i,j,k,iblock),
290 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
295 C For dummy ends assign glycine-type coefficients of theta-only terms; the
296 C coefficients of theta-and-gamma-dependent terms are zero.
297 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
298 C RECOMENTDED AFTER VERSION 3.3)
302 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
303 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
305 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
306 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
309 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
311 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
314 C AND COMMENT THE LOOPS BELOW
318 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
319 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
321 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
322 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
325 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
327 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
331 C Substitution for D aminoacids from symmetry.
334 do j=-nthetyp,nthetyp
335 do k=-nthetyp,nthetyp
336 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
338 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
342 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
343 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
344 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
345 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
351 ffthet(llll,lll,ll,i,j,k,iblock)=
352 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
353 ffthet(lll,llll,ll,i,j,k,iblock)=
354 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
355 ggthet(llll,lll,ll,i,j,k,iblock)=
356 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
357 ggthet(lll,llll,ll,i,j,k,iblock)=
358 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
367 C Control printout of the coefficients of virtual-bond-angle potentials
370 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
373 do j=-nthetyp,nthetyp
374 do k=-nthetyp,nthetyp
375 write (iout,'(//4a)')
376 & 'Type ',toronelet(i),toronelet(j),toronelet(k)
377 write (iout,'(//a,10x,a)') " l","a[l]"
378 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
379 write (iout,'(i2,1pe15.5)')
380 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
382 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
383 & "b",l,"c",l,"d",l,"e",l
385 write (iout,'(i2,4(1pe15.5))') m,
386 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
387 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
391 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
392 & "f+",l,"f-",l,"g+",l,"g-",l
395 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
396 & ffthet(n,m,l,i,j,k,iblock),
397 & ffthet(m,n,l,i,j,k,iblock),
398 & ggthet(n,m,l,i,j,k,iblock),
399 & ggthet(m,n,l,i,j,k,iblock)
412 C here will be the apropriate recalibrating for D-aminoacid
413 read (ithep,*,end=121,err=121) nthetyp
414 do i=-nthetyp+1,nthetyp-1
415 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
416 do j=0,nbend_kcc_Tb(i)
417 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
422 & "Parameters of the valence-only potentials"
423 do i=-nthetyp+1,nthetyp-1
424 write (iout,'(2a)') "Type ",toronelet(i)
425 do k=0,nbend_kcc_Tb(i)
426 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
433 c write (2,*) "Start reading THETA_PDB",ithep_pdb
436 read (ithep_pdb,*,err=111,end=111)
437 & a0thet(i),(athet(j,i,1,1),j=1,2),
438 & (bthet(j,i,1,1),j=1,2)
439 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
440 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
441 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
445 athet(1,i,1,-1)=athet(1,i,1,1)
446 athet(2,i,1,-1)=athet(2,i,1,1)
447 bthet(1,i,1,-1)=-bthet(1,i,1,1)
448 bthet(2,i,1,-1)=-bthet(2,i,1,1)
449 athet(1,i,-1,1)=-athet(1,i,1,1)
450 athet(2,i,-1,1)=-athet(2,i,1,1)
451 bthet(1,i,-1,1)=bthet(1,i,1,1)
452 bthet(2,i,-1,1)=bthet(2,i,1,1)
456 athet(1,i,-1,-1)=athet(1,-i,1,1)
457 athet(2,i,-1,-1)=-athet(2,-i,1,1)
458 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
459 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
460 athet(1,i,-1,1)=athet(1,-i,1,1)
461 athet(2,i,-1,1)=-athet(2,-i,1,1)
462 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
463 bthet(2,i,-1,1)=bthet(2,-i,1,1)
464 athet(1,i,1,-1)=-athet(1,-i,1,1)
465 athet(2,i,1,-1)=athet(2,-i,1,1)
466 bthet(1,i,1,-1)=bthet(1,-i,1,1)
467 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
472 polthet(j,i)=polthet(j,-i)
475 gthet(j,i)=gthet(j,-i)
478 c write (2,*) "End reading THETA_PDB"
484 C Read the parameters of the probability distribution/energy expression
485 C of the side chains.
488 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
492 dsc_inv(i)=1.0D0/dsc(i)
503 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
504 & ((blower(k,l,1),l=1,k),k=1,3)
505 censc(1,1,-i)=censc(1,1,i)
506 censc(2,1,-i)=censc(2,1,i)
507 censc(3,1,-i)=-censc(3,1,i)
509 read (irotam,*,end=112,err=112) bsc(j,i)
510 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
511 & ((blower(k,l,j),l=1,k),k=1,3)
512 censc(1,j,-i)=censc(1,j,i)
513 censc(2,j,-i)=censc(2,j,i)
514 censc(3,j,-i)=-censc(3,j,i)
515 C BSC is amplitude of Gaussian
522 akl=akl+blower(k,m,j)*blower(l,m,j)
526 if (((k.eq.3).and.(l.ne.3))
527 & .or.((l.eq.3).and.(k.ne.3))) then
528 gaussc(k,l,j,-i)=-akl
529 gaussc(l,k,j,-i)=-akl
541 write (iout,'(/a)') 'Parameters of side-chain local geometry'
546 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
547 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
548 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
549 & 'log h',(bsc(j,i),j=1,nlobi)
550 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
551 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
553 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
554 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
557 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
558 write (iout,'(a,f10.4,4(16x,f10.4))')
559 & 'Center ',(bsc(j,i),j=1,nlobi)
560 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
569 C Read scrot parameters for potentials determined from all-atom AM1 calculations
570 C added by Urszula Kozlowska 07/11/2007
573 read (irotam,*,end=112,err=112)
575 read (irotam,*,end=112,err=112)
578 read(irotam,*,end=112,err=112) sc_parmin(j,i)
583 C Read the parameters of the probability distribution/energy expression
584 C of the side chains.
586 write (2,*) "Start reading ROTAM_PDB"
588 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
592 dsc_inv(i)=1.0D0/dsc(i)
603 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
604 & ((blower(k,l,1),l=1,k),k=1,3)
606 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
607 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
608 & ((blower(k,l,j),l=1,k),k=1,3)
615 akl=akl+blower(k,m,j)*blower(l,m,j)
625 write (2,*) "End reading ROTAM_PDB"
629 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
630 C interaction energy of the Gly, Ala, and Pro prototypes.
632 read (ifourier,*) nloctyp
633 SPLIT_FOURIERTOR = nloctyp.lt.0
634 nloctyp = iabs(nloctyp)
636 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
637 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
638 itype2loc(ntyp1)=nloctyp
639 iloctyp(nloctyp)=ntyp1
641 itype2loc(-i)=-itype2loc(i)
650 iloctyp(-i)=-iloctyp(i)
652 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
653 c write (iout,*) "nloctyp",nloctyp,
654 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
657 c write (iout,*) "NEWCORR",i
658 read (ifourier,*,end=115,err=115)
661 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
664 c write (iout,*) "NEWCORR BNEW1"
665 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
668 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
671 c write (iout,*) "NEWCORR BNEW2"
672 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
674 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
675 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
677 c write (iout,*) "NEWCORR CCNEW"
678 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
680 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
681 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
683 c write (iout,*) "NEWCORR DDNEW"
684 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
688 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
692 c write (iout,*) "NEWCORR EENEW1"
693 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
695 read (ifourier,*,end=115,err=115) e0new(ii,i)
697 c write (iout,*) (e0new(ii,i),ii=1,3)
699 c write (iout,*) "NEWCORR EENEW"
702 ccnew(ii,1,i)=ccnew(ii,1,i)/2
703 ccnew(ii,2,i)=ccnew(ii,2,i)/2
704 ddnew(ii,1,i)=ddnew(ii,1,i)/2
705 ddnew(ii,2,i)=ddnew(ii,2,i)/2
710 bnew1(ii,1,-i)= bnew1(ii,1,i)
711 bnew1(ii,2,-i)=-bnew1(ii,2,i)
712 bnew2(ii,1,-i)= bnew2(ii,1,i)
713 bnew2(ii,2,-i)=-bnew2(ii,2,i)
716 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
717 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
718 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
719 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
720 ccnew(ii,1,-i)=ccnew(ii,1,i)
721 ccnew(ii,2,-i)=-ccnew(ii,2,i)
722 ddnew(ii,1,-i)=ddnew(ii,1,i)
723 ddnew(ii,2,-i)=-ddnew(ii,2,i)
725 e0new(1,-i)= e0new(1,i)
726 e0new(2,-i)=-e0new(2,i)
727 e0new(3,-i)=-e0new(3,i)
729 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
730 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
731 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
732 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
736 write (iout,'(a)') "Coefficients of the multibody terms"
737 do i=-nloctyp+1,nloctyp-1
738 write (iout,*) "Type: ",onelet(iloctyp(i))
739 write (iout,*) "Coefficients of the expansion of B1"
741 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
743 write (iout,*) "Coefficients of the expansion of B2"
745 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
747 write (iout,*) "Coefficients of the expansion of C"
748 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
749 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
750 write (iout,*) "Coefficients of the expansion of D"
751 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
752 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
753 write (iout,*) "Coefficients of the expansion of E"
754 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
757 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
762 IF (SPLIT_FOURIERTOR) THEN
764 c write (iout,*) "NEWCORR TOR",i
765 read (ifourier,*,end=115,err=115)
768 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
771 c write (iout,*) "NEWCORR BNEW1 TOR"
772 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
775 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
778 c write (iout,*) "NEWCORR BNEW2 TOR"
779 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
781 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
782 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
784 c write (iout,*) "NEWCORR CCNEW TOR"
785 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
787 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
788 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
790 c write (iout,*) "NEWCORR DDNEW TOR"
791 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
795 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
799 c write (iout,*) "NEWCORR EENEW1 TOR"
800 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
802 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
804 c write (iout,*) (e0newtor(ii,i),ii=1,3)
806 c write (iout,*) "NEWCORR EENEW TOR"
809 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
810 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
811 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
812 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
817 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
818 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
819 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
820 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
823 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
824 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
825 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
826 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
828 e0newtor(1,-i)= e0newtor(1,i)
829 e0newtor(2,-i)=-e0newtor(2,i)
830 e0newtor(3,-i)=-e0newtor(3,i)
832 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
833 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
834 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
835 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
840 & "Single-body coefficients of the torsional potentials"
841 do i=-nloctyp+1,nloctyp-1
842 write (iout,*) "Type: ",onelet(iloctyp(i))
843 write (iout,*) "Coefficients of the expansion of B1tor"
845 write (iout,'(3hB1(,i1,1h),3f10.5)')
846 & j,(bnew1tor(k,j,i),k=1,3)
848 write (iout,*) "Coefficients of the expansion of B2tor"
850 write (iout,'(3hB2(,i1,1h),3f10.5)')
851 & j,(bnew2tor(k,j,i),k=1,3)
853 write (iout,*) "Coefficients of the expansion of Ctor"
854 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
855 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
856 write (iout,*) "Coefficients of the expansion of Dtor"
857 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
858 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
859 write (iout,*) "Coefficients of the expansion of Etor"
860 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
863 write (iout,'(1hE,2i1,2f10.5)')
864 & j,k,(eenewtor(l,j,k,i),l=1,2)
870 do i=-nloctyp+1,nloctyp-1
873 bnew1tor(ii,j,i)=bnew1(ii,j,i)
878 bnew2tor(ii,j,i)=bnew2(ii,j,i)
882 ccnewtor(ii,1,i)=ccnew(ii,1,i)
883 ccnewtor(ii,2,i)=ccnew(ii,2,i)
884 ddnewtor(ii,1,i)=ddnew(ii,1,i)
885 ddnewtor(ii,2,i)=ddnew(ii,2,i)
888 ENDIF !SPLIT_FOURIER_TOR
891 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
893 read (ifourier,*,end=115,err=115)
894 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
896 write (iout,*) 'Type ',onelet(iloctyp(i))
897 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
911 c B1tilde(1,i) = b(3)
912 c B1tilde(2,i) =-b(5)
913 c B1tilde(1,-i) =-b(3)
914 c B1tilde(2,-i) =b(5)
921 cc B1tilde(1,i) = b(3,i)
922 cc B1tilde(2,i) =-b(5,i)
923 C B1tilde(1,-i) =-b(3,i)
924 C B1tilde(2,-i) =b(5,i)
925 cc b1tilde(1,i)=0.0d0
926 cc b1tilde(2,i)=0.0d0
938 CCold(1,1,-i)= b(7,i)
939 CCold(2,2,-i)=-b(7,i)
940 CCold(2,1,-i)=-b(9,i)
941 CCold(1,2,-i)=-b(9,i)
946 c Ctilde(1,1,i)= CCold(1,1,i)
947 c Ctilde(1,2,i)= CCold(1,2,i)
948 c Ctilde(2,1,i)=-CCold(2,1,i)
949 c Ctilde(2,2,i)=-CCold(2,2,i)
951 c Ctilde(1,1,i)=0.0d0
952 c Ctilde(1,2,i)=0.0d0
953 c Ctilde(2,1,i)=0.0d0
954 c Ctilde(2,2,i)=0.0d0
959 DDold(1,1,-i)= b(6,i)
960 DDold(2,2,-i)=-b(6,i)
961 DDold(2,1,-i)=-b(8,i)
962 DDold(1,2,-i)=-b(8,i)
967 c Dtilde(1,1,i)= DD(1,1,i)
968 c Dtilde(1,2,i)= DD(1,2,i)
969 c Dtilde(2,1,i)=-DD(2,1,i)
970 c Dtilde(2,2,i)=-DD(2,2,i)
972 c Dtilde(1,1,i)=0.0d0
973 c Dtilde(1,2,i)=0.0d0
974 c Dtilde(2,1,i)=0.0d0
975 c Dtilde(2,2,i)=0.0d0
976 EEold(1,1,i)= b(10,i)+b(11,i)
977 EEold(2,2,i)=-b(10,i)+b(11,i)
978 EEold(2,1,i)= b(12,i)-b(13,i)
979 EEold(1,2,i)= b(12,i)+b(13,i)
980 EEold(1,1,-i)= b(10,i)+b(11,i)
981 EEold(2,2,-i)=-b(10,i)+b(11,i)
982 EEold(2,1,-i)=-b(12,i)+b(13,i)
983 EEold(1,2,-i)=-b(12,i)-b(13,i)
984 c write(iout,*) "TU DOCHODZE"
990 c ee(2,1,i)=ee(1,2,i)
995 &"Coefficients of the cumulants (independent of valence angles)"
996 do i=-nloctyp+1,nloctyp-1
997 write (iout,*) 'Type ',onelet(iloctyp(i))
999 write(iout,'(2f10.5)') B(3,i),B(5,i)
1001 write(iout,'(2f10.5)') B(2,i),B(4,i)
1004 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1008 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1012 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1020 C Read torsional parameters in old format
1022 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1023 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1024 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1029 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1035 write (iout,'(/a/)') 'Torsional constants:'
1038 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1039 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1045 C Read torsional parameters
1047 IF (TOR_MODE.eq.0) THEN
1049 read (itorp,*,end=113,err=113) ntortyp
1050 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1052 itype2loc(i)=itortyp(i)
1055 itype2loc(-i)=-itype2loc(i)
1057 itortyp(ntyp1)=ntortyp
1060 itortyp(i)=-itortyp(-i)
1062 write (iout,*) 'ntortyp',ntortyp
1064 do j=-ntortyp+1,ntortyp-1
1065 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1067 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1068 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1071 do k=1,nterm(i,j,iblock)
1072 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1074 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1075 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1076 v0ij=v0ij+si*v1(k,i,j,iblock)
1078 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
1079 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
1080 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
1082 do k=1,nlor(i,j,iblock)
1083 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1084 & vlor2(k,i,j),vlor3(k,i,j)
1085 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1088 v0(-i,-j,iblock)=v0ij
1094 write (iout,'(/a/)') 'Parameters of the SCCOR potentials:'
1097 do j=-ntortyp+1,ntortyp-1
1098 write (iout,*) 'ityp',i,' jtyp',j
1099 write (iout,*) 'Fourier constants'
1100 do k=1,nterm(i,j,iblock)
1101 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1104 write (iout,*) 'Lorenz constants'
1105 do k=1,nlor(i,j,iblock)
1106 write (iout,'(3(1pe15.5))')
1107 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1115 C 6/23/01 Read parameters for double torsionals
1119 do j=-ntortyp+1,ntortyp-1
1120 do k=-ntortyp+1,ntortyp-1
1121 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1122 c write (iout,*) "OK onelett",
1125 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1126 & .or. t3.ne.toronelet(k)) then
1127 write (iout,*) "Error in double torsional parameter file",
1130 call MPI_Finalize(Ierror)
1132 stop "Error in double torsional parameter file"
1134 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1135 & ntermd_2(i,j,k,iblock)
1136 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1137 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1138 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1139 & ntermd_1(i,j,k,iblock))
1140 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1141 & ntermd_1(i,j,k,iblock))
1142 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1143 & ntermd_1(i,j,k,iblock))
1144 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1145 & ntermd_1(i,j,k,iblock))
1146 C Martix of D parameters for one dimesional foureir series
1147 do l=1,ntermd_1(i,j,k,iblock)
1148 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1149 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1150 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1151 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1152 c write(iout,*) "whcodze" ,
1153 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1155 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1156 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1157 & v2s(m,l,i,j,k,iblock),
1158 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1159 C Martix of D parameters for two dimesional fourier series
1160 do l=1,ntermd_2(i,j,k,iblock)
1162 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1163 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1164 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1165 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1174 write (iout,*) 'Constants for double torsionals'
1177 do j=-ntortyp+1,ntortyp-1
1178 do k=-ntortyp+1,ntortyp-1
1179 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1180 & ' nsingle',ntermd_1(i,j,k,iblock),
1181 & ' ndouble',ntermd_2(i,j,k,iblock)
1183 write (iout,*) 'Single angles:'
1184 do l=1,ntermd_1(i,j,k,iblock)
1185 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1186 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1187 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1188 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1191 write (iout,*) 'Pairs of angles:'
1192 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1193 do l=1,ntermd_2(i,j,k,iblock)
1194 write (iout,'(i5,20f10.5)')
1195 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1198 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1199 do l=1,ntermd_2(i,j,k,iblock)
1200 write (iout,'(i5,20f10.5)')
1201 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1202 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1211 ELSE IF (TOR_MODE.eq.1) THEN
1213 C read valence-torsional parameters
1214 read (itorp,*,end=121,err=121) ntortyp
1216 write (iout,*) "Valence-torsional parameters read in ntortyp",
1218 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
1219 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1222 itype2loc(i)=itortyp(i)
1226 itortyp(i)=-itortyp(-i)
1228 do i=-ntortyp+1,ntortyp-1
1229 do j=-ntortyp+1,ntortyp-1
1230 C first we read the cos and sin gamma parameters
1231 read (itorp,'(13x,a)',end=121,err=121) string
1232 write (iout,*) i,j,string
1233 read (itorp,*,end=121,err=121)
1234 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1235 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1236 do k=1,nterm_kcc(j,i)
1237 do l=1,nterm_kcc_Tb(j,i)
1238 do ll=1,nterm_kcc_Tb(j,i)
1239 read (itorp,*,end=121,err=121) ii,jj,kk,
1240 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1248 c AL 4/8/16: Calculate coefficients from one-body parameters
1251 itortyp(i)=itype2loc(i)
1254 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1256 do i=-ntortyp+1,ntortyp-1
1257 do j=-ntortyp+1,ntortyp-1
1260 do k=1,nterm_kcc_Tb(j,i)
1261 do l=1,nterm_kcc_Tb(j,i)
1262 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1263 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1264 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1265 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1268 do k=1,nterm_kcc_Tb(j,i)
1269 do l=1,nterm_kcc_Tb(j,i)
1271 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1272 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1273 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1274 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1276 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1277 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1278 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1279 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1283 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)
1287 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1292 if (tor_mode.gt.0 .and. lprint) then
1293 c Print valence-torsional parameters
1295 & "Parameters of the valence-torsional potentials"
1296 do i=-ntortyp+1,ntortyp-1
1297 do j=-ntortyp+1,ntortyp-1
1298 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1299 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1300 do k=1,nterm_kcc(j,i)
1301 do l=1,nterm_kcc_Tb(j,i)
1302 do ll=1,nterm_kcc_Tb(j,i)
1303 write (iout,'(3i5,2f15.4)')
1304 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1313 C Read of Side-chain backbone correlation parameters
1314 C Modified 11 May 2012 by Adasko
1317 read (isccor,*,end=119,err=119) nsccortyp
1319 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1321 isccortyp(i)=-isccortyp(-i)
1323 iscprol=isccortyp(20)
1324 c write (iout,*) 'ntortyp',ntortyp
1326 cc maxinter is maximum interaction sites
1330 read (isccor,*,end=119,err=119)
1331 &nterm_sccor(i,j),nlor_sccor(i,j)
1337 nterm_sccor(-i,j)=nterm_sccor(i,j)
1338 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1339 nterm_sccor(i,-j)=nterm_sccor(i,j)
1340 do k=1,nterm_sccor(i,j)
1341 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1343 if (j.eq.iscprol) then
1344 if (i.eq.isccortyp(10)) then
1345 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1346 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1348 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1349 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1350 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1351 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1352 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1353 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1354 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1355 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1358 if (i.eq.isccortyp(10)) then
1359 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1360 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1362 if (j.eq.isccortyp(10)) then
1363 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1364 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1366 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1367 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)
1375 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1376 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1377 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1378 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1381 do k=1,nlor_sccor(i,j)
1382 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1383 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1384 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1385 &(1+vlor3sccor(k,i,j)**2)
1387 v0sccor(l,i,j)=v0ijsccor
1388 v0sccor(l,-i,j)=v0ijsccor1
1389 v0sccor(l,i,-j)=v0ijsccor2
1390 v0sccor(l,-i,-j)=v0ijsccor3
1396 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1397 c write (iout,*) 'ntortyp',ntortyp
1399 cc maxinter is maximum interaction sites
1403 read (isccor,*,end=119,err=119)
1404 & nterm_sccor(i,j),nlor_sccor(i,j)
1408 do k=1,nterm_sccor(i,j)
1409 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1411 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1414 do k=1,nlor_sccor(i,j)
1415 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1416 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1417 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1418 &(1+vlor3sccor(k,i,j)**2)
1420 v0sccor(l,i,j)=v0ijsccor
1428 write (iout,'(/a/)') 'Torsional constants:'
1432 write (iout,*) 'ityp',i,' jtyp',j
1433 write (iout,*) 'Fourier constants'
1434 do k=1,nterm_sccor(i,j)
1435 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1437 write (iout,*) 'Lorenz constants'
1438 do k=1,nlor_sccor(i,j)
1439 write (iout,'(3(1pe15.5))')
1440 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1448 C Read electrostatic-interaction parameters
1452 write (iout,'(/a)') 'Electrostatic interaction constants:'
1453 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1454 & 'IT','JT','APP','BPP','AEL6','AEL3'
1456 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1457 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1458 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1459 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1464 app (i,j)=epp(i,j)*rri*rri
1465 bpp (i,j)=-2.0D0*epp(i,j)*rri
1466 ael6(i,j)=elpp6(i,j)*4.2D0**6
1467 ael3(i,j)=elpp3(i,j)*4.2D0**3
1469 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1470 & ael6(i,j),ael3(i,j)
1475 C Read side-chain interaction parameters.
1477 read (isidep,*,end=117,err=117) ipot,expon
1478 if (ipot.lt.1 .or. ipot.gt.5) then
1479 write (iout,'(2a)') 'Error while reading SC interaction',
1480 & 'potential file - unknown potential type.'
1482 call MPI_Finalize(Ierror)
1488 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1489 & ', exponents are ',expon,2*expon
1490 goto (10,20,30,30,40) ipot
1491 C----------------------- LJ potential ---------------------------------
1492 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1493 & (sigma0(i),i=1,ntyp)
1495 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1496 write (iout,'(a/)') 'The epsilon array:'
1497 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1498 write (iout,'(/a)') 'One-body parameters:'
1499 write (iout,'(a,4x,a)') 'residue','sigma'
1500 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1503 C----------------------- LJK potential --------------------------------
1504 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1505 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1507 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1508 write (iout,'(a/)') 'The epsilon array:'
1509 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1510 write (iout,'(/a)') 'One-body parameters:'
1511 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1512 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1516 C---------------------- GB or BP potential -----------------------------
1518 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1520 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1521 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1522 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1523 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1524 C now we start reading lipid
1526 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1527 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1529 C print *,"WARNING!!"
1531 C epslip(i,j)=epslip(i,j)+0.05d0
1534 if (lprint) write(iout,*) epslip(1,1),"OK?"
1535 C For the GB potential convert sigma'**2 into chi'
1538 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1542 write (iout,'(/a/)') 'Parameters of the BP potential:'
1543 write (iout,'(a/)') 'The epsilon array:'
1544 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1545 write (iout,'(/a)') 'One-body parameters:'
1546 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1548 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1549 & chip(i),alp(i),i=1,ntyp)
1552 C--------------------- GBV potential -----------------------------------
1553 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1554 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1555 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1557 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1558 write (iout,'(a/)') 'The epsilon array:'
1559 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1560 write (iout,'(/a)') 'One-body parameters:'
1561 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1562 & 's||/s_|_^2',' chip ',' alph '
1563 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1564 & sigii(i),chip(i),alp(i),i=1,ntyp)
1568 C-----------------------------------------------------------------------
1569 C Calculate the "working" parameters of SC interactions.
1573 epslip(i,j)=epslip(j,i)
1578 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1579 sigma(j,i)=sigma(i,j)
1580 rs0(i,j)=dwa16*sigma(i,j)
1584 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1585 & 'Working parameters of the SC interactions:',
1586 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1591 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1600 sigeps=dsign(1.0D0,epsij)
1602 aa_aq(i,j)=epsij*rrij*rrij
1603 bb_aq(i,j)=-sigeps*epsij*rrij
1604 aa_aq(j,i)=aa_aq(i,j)
1605 bb_aq(j,i)=bb_aq(i,j)
1606 epsijlip=epslip(i,j)
1607 sigeps=dsign(1.0D0,epsijlip)
1608 epsijlip=dabs(epsijlip)
1609 aa_lip(i,j)=epsijlip*rrij*rrij
1610 bb_lip(i,j)=-sigeps*epsijlip*rrij
1611 aa_lip(j,i)=aa_lip(i,j)
1612 bb_lip(j,i)=bb_lip(i,j)
1614 sigt1sq=sigma0(i)**2
1615 sigt2sq=sigma0(j)**2
1618 ratsig1=sigt2sq/sigt1sq
1619 ratsig2=1.0D0/ratsig1
1620 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1621 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1622 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1626 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1627 sigmaii(i,j)=rsum_max
1628 sigmaii(j,i)=rsum_max
1630 c sigmaii(i,j)=r0(i,j)
1631 c sigmaii(j,i)=r0(i,j)
1633 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1634 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1635 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1636 augm(i,j)=epsij*r_augm**(2*expon)
1637 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1644 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1645 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1646 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1651 write(iout,*) "tube param"
1652 read(itube,*) epspeptube,sigmapeptube
1653 sigmapeptube=sigmapeptube**6
1654 sigeps=dsign(1.0D0,epspeptube)
1655 epspeptube=dabs(epspeptube)
1656 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1657 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1658 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1660 read(itube,*) epssctube,sigmasctube
1661 sigmasctube=sigmasctube**6
1662 sigeps=dsign(1.0D0,epssctube)
1663 epssctube=dabs(epssctube)
1664 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1665 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1666 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1672 C Define the SC-p interaction constants (hard-coded; old style)
1675 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1677 c aad(i,1)=0.3D0*4.0D0**12
1678 C Following line for constants currently implemented
1679 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1680 aad(i,1)=1.5D0*4.0D0**12
1681 c aad(i,1)=0.17D0*5.6D0**12
1683 C "Soft" SC-p repulsion
1685 C Following line for constants currently implemented
1686 c aad(i,1)=0.3D0*4.0D0**6
1687 C "Hard" SC-p repulsion
1688 bad(i,1)=3.0D0*4.0D0**6
1689 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1698 C 8/9/01 Read the SC-p interaction constants from file
1701 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1704 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1705 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1706 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1707 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1711 write (iout,*) "Parameters of SC-p interactions:"
1713 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1714 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1720 C Define the constants of the disulfide bridge
1724 c Old arbitrary potential - commented out.
1729 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1730 c energy surface of diethyl disulfide.
1731 c A. Liwo and U. Kozlowska, 11/24/03
1747 C if(me.eq.king) then
1748 C write (iout,'(/a)') "Disulfide bridge parameters:"
1749 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1750 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1751 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1752 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1755 if (shield_mode.gt.0) then
1756 C VSolvSphere the volume of solving sphere
1757 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1758 C there will be no distinction between proline peptide group and normal peptide
1759 C group in case of shielding parameters
1760 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
1761 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1762 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1763 write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
1765 C long axis of side chain
1767 long_r_sidechain(i)=vbldsc0(1,i)
1768 short_r_sidechain(i)=sigma0(i)
1773 111 write (iout,*) "Error reading bending energy parameters."
1775 112 write (iout,*) "Error reading rotamer energy parameters."
1777 113 write (iout,*) "Error reading torsional energy parameters."
1779 114 write (iout,*) "Error reading double torsional energy parameters."
1782 & "Error reading cumulant (multibody energy) parameters."
1784 116 write (iout,*) "Error reading electrostatic energy parameters."
1786 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1788 117 write (iout,*) "Error reading side chain interaction parameters."
1790 118 write (iout,*) "Error reading SCp interaction parameters."
1792 119 write (iout,*) "Error reading SCCOR parameters"
1794 121 write (iout,*) "Error in Czybyshev parameters"
1797 call MPI_Finalize(Ierror)
1804 subroutine getenv_loc(var, val)
1805 character(*) var, val
1808 character(2000) line
1811 open (196,file='env',status='old',readonly,shared)
1813 c write(*,*)'looking for ',var
1814 10 read(196,*,err=11,end=11)line
1815 iread=index(line,var)
1816 c write(*,*)iread,' ',var,' ',line
1817 if (iread.eq.0) go to 10
1818 c write(*,*)'---> ',line
1824 iread=iread+ilen(var)+1
1825 read (line(iread:),*,err=12,end=12) val
1826 c write(*,*)'OK: ',var,' = ',val
1832 #elif (defined CRAY)
1833 integer lennam,lenval,ierror
1835 c getenv using a POSIX call, useful on the T3D
1836 c Sept 1996, comment out error check on advice of H. Pritchard
1839 if(lennam.le.0) stop '--error calling getenv--'
1840 call pxfgetenv(var,lennam,val,lenval,ierror)
1841 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1843 call getenv(var,val)
1845 C set the variables used for shielding effect
1846 C if (shield_mode.gt.0) then
1847 C VSolvSphere the volume of solving sphere
1849 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1850 C there will be no distinction between proline peptide group and normal peptide
1851 C group in case of shielding parameters
1852 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1853 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1854 C long axis of side chain
1856 C long_r_sidechain(i)=vbldsc0(1,i)
1857 C short_r_sidechain(i)=sigma0(i)
1859 C lets set the buffor value
1864 c---------------------------------------------------------------------------
1865 subroutine weightread
1867 C Read the parameters of the probability distributions of the virtual-bond
1868 C valence angles and the side chains and energy parameters.
1870 C Important! Energy-term weights ARE NOT read here; they are read from the
1871 C main input file instead, because NO defaults have yet been set for these
1874 implicit real*8 (a-h,o-z)
1875 include 'DIMENSIONS'
1880 include 'COMMON.IOUNITS'
1881 include 'COMMON.CHAIN'
1882 include 'COMMON.INTERACT'
1883 include 'COMMON.GEO'
1884 include 'COMMON.LOCAL'
1885 include 'COMMON.TORSION'
1886 include 'COMMON.SCCOR'
1887 include 'COMMON.SCROT'
1888 include 'COMMON.FFIELD'
1889 include 'COMMON.NAMES'
1890 include 'COMMON.SBRIDGE'
1892 include 'COMMON.SETUP'
1893 include 'COMMON.CONTROL'
1894 include 'COMMON.SHIELD'
1895 character*1000 weightcard
1897 c READ energy-term weights
1899 call card_concat(weightcard)
1900 call reada(weightcard,'WLONG',wlong,1.0D0)
1901 call reada(weightcard,'WSC',wsc,wlong)
1902 call reada(weightcard,'WSCP',wscp,wlong)
1903 call reada(weightcard,'WELEC',welec,1.0D0)
1904 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1905 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1906 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1907 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1908 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1909 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1910 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1911 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1912 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1913 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1914 call reada(weightcard,'WBOND',wbond,1.0D0)
1915 call reada(weightcard,'WTOR',wtor,1.0D0)
1916 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1917 call reada(weightcard,'WANG',wang,1.0D0)
1918 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1919 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
1920 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
1921 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
1922 call reada(weightcard,'WDFAB',wdfa_beta,0.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
1956 weights(28)=wdfa_dist
1957 weights(29)=wdfa_tor
1958 weights(30)=wdfa_nei
1959 weights(31)=wdfa_beta
1961 aad(i,1)=scalscp*aad(i,1)
1962 aad(i,2)=scalscp*aad(i,2)
1963 bad(i,1)=scalscp*bad(i,1)
1964 bad(i,2)=scalscp*bad(i,2)
1966 if(me.eq.king.or..not.out1file)
1967 & write (iout,100) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
1968 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
1970 100 format (/'Energy-term weights (unscaled):'//
1971 & 'WSCC= ',f10.6,' (SC-SC)'/
1972 & 'WSCP= ',f10.6,' (SC-p)'/
1973 & 'WELEC= ',f10.6,' (p-p electr)'/
1974 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
1975 & 'WBOND= ',f10.6,' (stretching)'/
1976 & 'WANG= ',f10.6,' (bending)'/
1977 & 'WSCLOC= ',f10.6,' (SC local)'/
1978 & 'WTOR= ',f10.6,' (torsional)'/
1979 & 'WTORD= ',f10.6,' (double torsional)'/
1980 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
1981 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
1982 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
1983 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
1984 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
1985 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
1986 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
1987 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
1988 & 'WTURN6= ',f10.6,' (turns, 6th order)')
1989 if(me.eq.king.or..not.out1file)then
1990 if (wcorr4.gt.0.0d0) then
1991 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
1992 & 'between contact pairs of peptide groups'
1993 write (iout,'(2(a,f5.3/))')
1994 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
1995 & 'Range of quenching the correlation terms:',2*delt_corr
1996 else if (wcorr.gt.0.0d0) then
1997 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
1998 & 'between contact pairs of peptide groups'
2000 write (iout,'(a,f8.3)')
2001 & 'Scaling factor of 1,4 SC-p interactions:',scal14
2002 write (iout,'(a,f8.3)')
2003 & 'General scaling factor of SC-p interactions:',scalscp
2005 r0_corr=cutoff_corr-delt_corr
2007 call rescale_weights(t_bath)
2008 if(me.eq.king.or..not.out1file)
2009 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
2010 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
2012 22 format (/'Energy-term weights (scaled):'//
2013 & 'WSCC= ',f10.6,' (SC-SC)'/
2014 & 'WSCP= ',f10.6,' (SC-p)'/
2015 & 'WELEC= ',f10.6,' (p-p electr)'/
2016 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
2017 & 'WBOND= ',f10.6,' (stretching)'/
2018 & 'WANG= ',f10.6,' (bending)'/
2019 & 'WSCLOC= ',f10.6,' (SC local)'/
2020 & 'WTOR= ',f10.6,' (torsional)'/
2021 & 'WTORD= ',f10.6,' (double torsional)'/
2022 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
2023 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
2024 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
2025 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
2026 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
2027 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
2028 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
2029 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
2030 & 'WTURN6= ',f10.6,' (turns, 6th order)')
2031 if(me.eq.king.or..not.out1file)
2032 & write (iout,*) "Reference temperature for weights calculation:",
2035 write (iout,'(/a/)') "DFA pseudopotential parameters:"
2036 write (iout,'(a,f10.6,a)')
2037 & "WDFAD= ",wdfa_dist," (distance)",
2038 & "WDFAT= ",wdfa_tor," (backbone angles)",
2039 & "WDFAN= ",wdfa_nei," (neighbors)",
2040 & "WDFAB= ",wdfa_beta," (beta structure)"
2042 call reada(weightcard,"D0CM",d0cm,3.78d0)
2043 call reada(weightcard,"AKCM",akcm,15.1d0)
2044 call reada(weightcard,"AKTH",akth,11.0d0)
2045 call reada(weightcard,"AKCT",akct,12.0d0)
2046 call reada(weightcard,"V1SS",v1ss,-1.08d0)
2047 call reada(weightcard,"V2SS",v2ss,7.61d0)
2048 call reada(weightcard,"V3SS",v3ss,13.7d0)
2049 call reada(weightcard,"EBR",ebr,-5.50D0)
2050 call reada(weightcard,"ATRISS",atriss,0.301D0)
2051 call reada(weightcard,"BTRISS",btriss,0.021D0)
2052 call reada(weightcard,"CTRISS",ctriss,1.001D0)
2053 call reada(weightcard,"DTRISS",dtriss,1.001D0)
2054 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
2056 dyn_ss_mask(i)=.false.
2060 dyn_ssbond_ij(i,j)=1.0d300
2063 call reada(weightcard,"HT",Ht,0.0D0)
2065 ss_depth=ebr/wsc-0.25*eps(1,1)
2066 Ht=Ht/wsc-0.25*eps(1,1)
2067 akcm=akcm*wstrain/wsc
2068 akth=akth*wstrain/wsc
2069 akct=akct*wstrain/wsc
2070 v1ss=v1ss*wstrain/wsc
2071 v2ss=v2ss*wstrain/wsc
2072 v3ss=v3ss*wstrain/wsc
2074 if (wstrain.ne.0.0) then
2075 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
2080 write (iout,*) "wshield,", wshield
2081 if(me.eq.king.or..not.out1file) then
2082 write (iout,*) "Parameters of the SS-bond potential:"
2083 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
2085 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
2086 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
2087 write (iout,*)" HT",Ht
2088 write (iout,*) "Parameters of the 'trisulfide' potential"
2089 write (iout,*) "ATRISS=", atriss
2090 write (iout,*) "BTRISS=", btriss
2091 write (iout,*) "CTRISS=", ctriss
2092 write (iout,*) "DTRISS=", dtriss
2093 print *,'indpdb=',indpdb,' pdbref=',pdbref