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
40 C For printing parameters after they are read set the following in the UNRES
43 C setenv PRINT_PARM YES
45 C To print parameters in LaTeX format rather than as ASCII tables:
49 call getenv_loc("PRINT_PARM",lancuch)
50 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
51 call getenv_loc("LATEX",lancuch)
52 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
54 dwa16=2.0d0**(1.0d0/6.0d0)
56 C Assign virtual-bond length
61 c Read the virtual-bond parameters, masses, and moments of inertia
62 c and Stokes radii of the peptide group and side chains
65 read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
68 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
73 dsc_inv(i)=1.0D0/dsc(i)
77 read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
79 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
80 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
85 dsc_inv(i)=1.0D0/dsc(i)
90 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
91 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
93 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
95 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
96 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
98 write (iout,'(13x,3f10.5)')
99 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
103 C reading lipid parameters
105 write (iout,*) "iliptranpar",iliptranpar
108 read(iliptranpar,*) pepliptran
110 read(iliptranpar,*) liptranene(i)
115 C Read the parameters of the probability distribution/energy expression
116 C of the virtual-bond valence angles theta
119 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
120 & (bthet(j,i,1,1),j=1,2)
121 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
122 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
123 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
127 athet(1,i,1,-1)=athet(1,i,1,1)
128 athet(2,i,1,-1)=athet(2,i,1,1)
129 bthet(1,i,1,-1)=-bthet(1,i,1,1)
130 bthet(2,i,1,-1)=-bthet(2,i,1,1)
131 athet(1,i,-1,1)=-athet(1,i,1,1)
132 athet(2,i,-1,1)=-athet(2,i,1,1)
133 bthet(1,i,-1,1)=bthet(1,i,1,1)
134 bthet(2,i,-1,1)=bthet(2,i,1,1)
138 athet(1,i,-1,-1)=athet(1,-i,1,1)
139 athet(2,i,-1,-1)=-athet(2,-i,1,1)
140 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
141 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
142 athet(1,i,-1,1)=athet(1,-i,1,1)
143 athet(2,i,-1,1)=-athet(2,-i,1,1)
144 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
145 bthet(2,i,-1,1)=bthet(2,-i,1,1)
146 athet(1,i,1,-1)=-athet(1,-i,1,1)
147 athet(2,i,1,-1)=athet(2,-i,1,1)
148 bthet(1,i,1,-1)=bthet(1,-i,1,1)
149 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
154 polthet(j,i)=polthet(j,-i)
157 gthet(j,i)=gthet(j,-i)
165 & 'Parameters of the virtual-bond valence angles:'
166 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
167 & ' ATHETA0 ',' A1 ',' A2 ',
170 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
171 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
173 write (iout,'(/a/9x,5a/79(1h-))')
174 & 'Parameters of the expression for sigma(theta_c):',
175 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
176 & ' ALPH3 ',' SIGMA0C '
178 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
179 & (polthet(j,i),j=0,3),sigc0(i)
181 write (iout,'(/a/9x,5a/79(1h-))')
182 & 'Parameters of the second gaussian:',
183 & ' THETA0 ',' SIGMA0 ',' G1 ',
186 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
187 & sig0(i),(gthet(j,i),j=1,3)
191 & 'Parameters of the virtual-bond valence angles:'
192 write (iout,'(/a/9x,5a/79(1h-))')
193 & 'Coefficients of expansion',
194 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
195 & ' b1*10^1 ',' b2*10^1 '
197 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
198 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
199 & (10*bthet(j,i,1,1),j=1,2)
201 write (iout,'(/a/9x,5a/79(1h-))')
202 & 'Parameters of the expression for sigma(theta_c):',
203 & ' alpha0 ',' alph1 ',' alph2 ',
204 & ' alhp3 ',' sigma0c '
206 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
207 & (polthet(j,i),j=0,3),sigc0(i)
209 write (iout,'(/a/9x,5a/79(1h-))')
210 & 'Parameters of the second gaussian:',
211 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
214 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
215 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
220 IF (tor_mode.eq.0) THEN
222 C Read the parameters of Utheta determined from ab initio surfaces
223 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
225 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
226 & ntheterm3,nsingle,ndouble
227 write (iout,*) "ithep",ithep
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
634 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
635 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
636 itype2loc(ntyp1)=nloctyp
637 iloctyp(nloctyp)=ntyp1
639 itype2loc(-i)=-itype2loc(i)
648 iloctyp(-i)=-iloctyp(i)
650 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
651 c write (iout,*) "nloctyp",nloctyp,
652 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
655 c write (iout,*) "NEWCORR",i
656 read (ifourier,*,end=115,err=115)
659 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
662 c write (iout,*) "NEWCORR BNEW1"
663 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
666 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
669 c write (iout,*) "NEWCORR BNEW2"
670 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
672 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
673 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
675 c write (iout,*) "NEWCORR CCNEW"
676 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
678 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
679 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
681 c write (iout,*) "NEWCORR DDNEW"
682 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
686 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
690 c write (iout,*) "NEWCORR EENEW1"
691 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
693 read (ifourier,*,end=115,err=115) e0new(ii,i)
695 c write (iout,*) (e0new(ii,i),ii=1,3)
697 c write (iout,*) "NEWCORR EENEW"
700 ccnew(ii,1,i)=ccnew(ii,1,i)/2
701 ccnew(ii,2,i)=ccnew(ii,2,i)/2
702 ddnew(ii,1,i)=ddnew(ii,1,i)/2
703 ddnew(ii,2,i)=ddnew(ii,2,i)/2
708 bnew1(ii,1,-i)= bnew1(ii,1,i)
709 bnew1(ii,2,-i)=-bnew1(ii,2,i)
710 bnew2(ii,1,-i)= bnew2(ii,1,i)
711 bnew2(ii,2,-i)=-bnew2(ii,2,i)
714 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
715 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
716 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
717 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
718 ccnew(ii,1,-i)=ccnew(ii,1,i)
719 ccnew(ii,2,-i)=-ccnew(ii,2,i)
720 ddnew(ii,1,-i)=ddnew(ii,1,i)
721 ddnew(ii,2,-i)=-ddnew(ii,2,i)
723 e0new(1,-i)= e0new(1,i)
724 e0new(2,-i)=-e0new(2,i)
725 e0new(3,-i)=-e0new(3,i)
727 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
728 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
729 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
730 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
734 write (iout,'(a)') "Coefficients of the multibody terms"
735 do i=-nloctyp+1,nloctyp-1
736 write (iout,*) "Type: ",onelet(iloctyp(i))
737 write (iout,*) "Coefficients of the expansion of B1"
739 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
741 write (iout,*) "Coefficients of the expansion of B2"
743 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
745 write (iout,*) "Coefficients of the expansion of C"
746 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
747 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
748 write (iout,*) "Coefficients of the expansion of D"
749 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
750 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
751 write (iout,*) "Coefficients of the expansion of E"
752 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
755 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
762 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
764 read (ifourier,*,end=115,err=115)
765 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
767 write (iout,*) 'Type ',onelet(iloctyp(i))
768 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
782 c B1tilde(1,i) = b(3)
783 c B1tilde(2,i) =-b(5)
784 c B1tilde(1,-i) =-b(3)
785 c B1tilde(2,-i) =b(5)
792 cc B1tilde(1,i) = b(3,i)
793 cc B1tilde(2,i) =-b(5,i)
794 C B1tilde(1,-i) =-b(3,i)
795 C B1tilde(2,-i) =b(5,i)
796 cc b1tilde(1,i)=0.0d0
797 cc b1tilde(2,i)=0.0d0
809 CCold(1,1,-i)= b(7,i)
810 CCold(2,2,-i)=-b(7,i)
811 CCold(2,1,-i)=-b(9,i)
812 CCold(1,2,-i)=-b(9,i)
817 c Ctilde(1,1,i)= CCold(1,1,i)
818 c Ctilde(1,2,i)= CCold(1,2,i)
819 c Ctilde(2,1,i)=-CCold(2,1,i)
820 c Ctilde(2,2,i)=-CCold(2,2,i)
822 c Ctilde(1,1,i)=0.0d0
823 c Ctilde(1,2,i)=0.0d0
824 c Ctilde(2,1,i)=0.0d0
825 c Ctilde(2,2,i)=0.0d0
830 DDold(1,1,-i)= b(6,i)
831 DDold(2,2,-i)=-b(6,i)
832 DDold(2,1,-i)=-b(8,i)
833 DDold(1,2,-i)=-b(8,i)
838 c Dtilde(1,1,i)= DD(1,1,i)
839 c Dtilde(1,2,i)= DD(1,2,i)
840 c Dtilde(2,1,i)=-DD(2,1,i)
841 c Dtilde(2,2,i)=-DD(2,2,i)
843 c Dtilde(1,1,i)=0.0d0
844 c Dtilde(1,2,i)=0.0d0
845 c Dtilde(2,1,i)=0.0d0
846 c Dtilde(2,2,i)=0.0d0
847 EEold(1,1,i)= b(10,i)+b(11,i)
848 EEold(2,2,i)=-b(10,i)+b(11,i)
849 EEold(2,1,i)= b(12,i)-b(13,i)
850 EEold(1,2,i)= b(12,i)+b(13,i)
851 EEold(1,1,-i)= b(10,i)+b(11,i)
852 EEold(2,2,-i)=-b(10,i)+b(11,i)
853 EEold(2,1,-i)=-b(12,i)+b(13,i)
854 EEold(1,2,-i)=-b(12,i)-b(13,i)
855 write(iout,*) "TU DOCHODZE"
861 c ee(2,1,i)=ee(1,2,i)
866 &"Coefficients of the cumulants (independent of valence angles)"
867 do i=-nloctyp+1,nloctyp-1
868 write (iout,*) 'Type ',onelet(iloctyp(i))
870 write(iout,'(2f10.5)') B(3,i),B(5,i)
872 write(iout,'(2f10.5)') B(2,i),B(4,i)
875 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
879 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
883 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
891 C Read torsional parameters in old format
893 read (itorp,*,end=113,err=113) ntortyp,nterm_old
894 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
895 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
900 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
906 write (iout,'(/a/)') 'Torsional constants:'
909 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
910 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
916 C Read torsional parameters
918 IF (TOR_MODE.eq.0) THEN
920 read (itorp,*,end=113,err=113) ntortyp
921 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
923 itype2loc(i)=itortyp(i)
926 itype2loc(-i)=-itype2loc(i)
928 itortyp(ntyp1)=ntortyp
931 itortyp(i)=-itortyp(-i)
933 write (iout,*) 'ntortyp',ntortyp
935 do j=-ntortyp+1,ntortyp-1
936 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
938 nterm(-i,-j,iblock)=nterm(i,j,iblock)
939 nlor(-i,-j,iblock)=nlor(i,j,iblock)
942 do k=1,nterm(i,j,iblock)
943 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
945 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
946 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
947 v0ij=v0ij+si*v1(k,i,j,iblock)
949 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
950 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
951 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
953 do k=1,nlor(i,j,iblock)
954 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
955 & vlor2(k,i,j),vlor3(k,i,j)
956 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
959 v0(-i,-j,iblock)=v0ij
965 write (iout,'(/a/)') 'Parameters of the SCCOR potentials:'
968 do j=-ntortyp+1,ntortyp-1
969 write (iout,*) 'ityp',i,' jtyp',j
970 write (iout,*) 'Fourier constants'
971 do k=1,nterm(i,j,iblock)
972 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
975 write (iout,*) 'Lorenz constants'
976 do k=1,nlor(i,j,iblock)
977 write (iout,'(3(1pe15.5))')
978 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
986 C 6/23/01 Read parameters for double torsionals
990 do j=-ntortyp+1,ntortyp-1
991 do k=-ntortyp+1,ntortyp-1
992 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
993 c write (iout,*) "OK onelett",
996 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
997 & .or. t3.ne.toronelet(k)) then
998 write (iout,*) "Error in double torsional parameter file",
1001 call MPI_Finalize(Ierror)
1003 stop "Error in double torsional parameter file"
1005 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1006 & ntermd_2(i,j,k,iblock)
1007 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1008 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1009 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1010 & ntermd_1(i,j,k,iblock))
1011 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1012 & ntermd_1(i,j,k,iblock))
1013 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1014 & ntermd_1(i,j,k,iblock))
1015 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1016 & ntermd_1(i,j,k,iblock))
1017 C Martix of D parameters for one dimesional foureir series
1018 do l=1,ntermd_1(i,j,k,iblock)
1019 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1020 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1021 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1022 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1023 c write(iout,*) "whcodze" ,
1024 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1026 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1027 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1028 & v2s(m,l,i,j,k,iblock),
1029 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1030 C Martix of D parameters for two dimesional fourier series
1031 do l=1,ntermd_2(i,j,k,iblock)
1033 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1034 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1035 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1036 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1045 write (iout,*) 'Constants for double torsionals'
1048 do j=-ntortyp+1,ntortyp-1
1049 do k=-ntortyp+1,ntortyp-1
1050 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1051 & ' nsingle',ntermd_1(i,j,k,iblock),
1052 & ' ndouble',ntermd_2(i,j,k,iblock)
1054 write (iout,*) 'Single angles:'
1055 do l=1,ntermd_1(i,j,k,iblock)
1056 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1057 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1058 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1059 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1062 write (iout,*) 'Pairs of angles:'
1063 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1064 do l=1,ntermd_2(i,j,k,iblock)
1065 write (iout,'(i5,20f10.5)')
1066 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1069 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1070 do l=1,ntermd_2(i,j,k,iblock)
1071 write (iout,'(i5,20f10.5)')
1072 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1073 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1082 ELSE IF (TOR_MODE.eq.1) THEN
1084 C read valence-torsional parameters
1085 read (itorp,*,end=121,err=121) ntortyp
1087 write (iout,*) "Valence-torsional parameters read in ntortyp",
1089 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
1090 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1093 itype2loc(i)=itortyp(i)
1097 itortyp(i)=-itortyp(-i)
1099 do i=-ntortyp+1,ntortyp-1
1100 do j=-ntortyp+1,ntortyp-1
1101 C first we read the cos and sin gamma parameters
1102 read (itorp,'(13x,a)',end=121,err=121) string
1103 write (iout,*) i,j,string
1104 read (itorp,*,end=121,err=121)
1105 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1106 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1107 do k=1,nterm_kcc(j,i)
1108 do l=1,nterm_kcc_Tb(j,i)
1109 do ll=1,nterm_kcc_Tb(j,i)
1110 read (itorp,*,end=121,err=121) ii,jj,kk,
1111 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1119 c AL 4/8/16: Calculate coefficients from one-body parameters
1122 itortyp(i)=itype2loc(i)
1125 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1127 do i=-ntortyp+1,ntortyp-1
1128 do j=-ntortyp+1,ntortyp-1
1131 do k=1,nterm_kcc_Tb(j,i)
1132 do l=1,nterm_kcc_Tb(j,i)
1133 v1_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,1,j)
1134 & +bnew1(k,2,i)*bnew2(l,2,j)
1135 v2_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,2,j)
1136 & +bnew1(k,2,i)*bnew2(l,1,j)
1139 do k=1,nterm_kcc_Tb(j,i)
1140 do l=1,nterm_kcc_Tb(j,i)
1141 v1_kcc(k,l,2,i,j)=-0.25*(ccnew(k,1,i)*ddnew(l,1,j)
1142 & -ccnew(k,2,i)*ddnew(l,2,j))
1143 v2_kcc(k,l,2,i,j)=-0.25*(ccnew(k,2,i)*ddnew(l,1,j)
1144 & +ccnew(k,1,i)*ddnew(l,2,j))
1147 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)
1151 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1156 if (tor_mode.gt.0 .and. lprint) then
1157 c Print valence-torsional parameters
1159 & "Parameters of the valence-torsional potentials"
1160 do i=-ntortyp+1,ntortyp-1
1161 do j=-ntortyp+1,ntortyp-1
1162 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1163 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1164 do k=1,nterm_kcc(j,i)
1165 do l=1,nterm_kcc_Tb(j,i)
1166 do ll=1,nterm_kcc_Tb(j,i)
1167 write (iout,'(3i5,2f15.4)')
1168 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1177 C Read of Side-chain backbone correlation parameters
1178 C Modified 11 May 2012 by Adasko
1181 read (isccor,*,end=119,err=119) nsccortyp
1183 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1185 isccortyp(i)=-isccortyp(-i)
1187 iscprol=isccortyp(20)
1188 c write (iout,*) 'ntortyp',ntortyp
1190 cc maxinter is maximum interaction sites
1194 read (isccor,*,end=119,err=119)
1195 &nterm_sccor(i,j),nlor_sccor(i,j)
1201 nterm_sccor(-i,j)=nterm_sccor(i,j)
1202 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1203 nterm_sccor(i,-j)=nterm_sccor(i,j)
1204 do k=1,nterm_sccor(i,j)
1205 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1207 if (j.eq.iscprol) then
1208 if (i.eq.isccortyp(10)) then
1209 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1210 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1212 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1213 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1214 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1215 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1216 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1217 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1218 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1219 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1222 if (i.eq.isccortyp(10)) then
1223 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1224 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1226 if (j.eq.isccortyp(10)) then
1227 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1228 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1230 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1231 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1232 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1233 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1234 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1235 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1239 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1240 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1241 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1242 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1245 do k=1,nlor_sccor(i,j)
1246 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1247 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1248 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1249 &(1+vlor3sccor(k,i,j)**2)
1251 v0sccor(l,i,j)=v0ijsccor
1252 v0sccor(l,-i,j)=v0ijsccor1
1253 v0sccor(l,i,-j)=v0ijsccor2
1254 v0sccor(l,-i,-j)=v0ijsccor3
1260 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1261 c write (iout,*) 'ntortyp',ntortyp
1263 cc maxinter is maximum interaction sites
1267 read (isccor,*,end=119,err=119)
1268 & nterm_sccor(i,j),nlor_sccor(i,j)
1272 do k=1,nterm_sccor(i,j)
1273 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1275 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1278 do k=1,nlor_sccor(i,j)
1279 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1280 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1281 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1282 &(1+vlor3sccor(k,i,j)**2)
1284 v0sccor(l,i,j)=v0ijsccor
1292 write (iout,'(/a/)') 'Torsional constants:'
1296 write (iout,*) 'ityp',i,' jtyp',j
1297 write (iout,*) 'Fourier constants'
1298 do k=1,nterm_sccor(i,j)
1299 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1301 write (iout,*) 'Lorenz constants'
1302 do k=1,nlor_sccor(i,j)
1303 write (iout,'(3(1pe15.5))')
1304 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1312 C Read electrostatic-interaction parameters
1316 write (iout,'(/a)') 'Electrostatic interaction constants:'
1317 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1318 & 'IT','JT','APP','BPP','AEL6','AEL3'
1320 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1321 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1322 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1323 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1328 app (i,j)=epp(i,j)*rri*rri
1329 bpp (i,j)=-2.0D0*epp(i,j)*rri
1330 ael6(i,j)=elpp6(i,j)*4.2D0**6
1331 ael3(i,j)=elpp3(i,j)*4.2D0**3
1333 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1334 & ael6(i,j),ael3(i,j)
1339 C Read side-chain interaction parameters.
1341 read (isidep,*,end=117,err=117) ipot,expon
1342 if (ipot.lt.1 .or. ipot.gt.5) then
1343 write (iout,'(2a)') 'Error while reading SC interaction',
1344 & 'potential file - unknown potential type.'
1346 call MPI_Finalize(Ierror)
1352 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1353 & ', exponents are ',expon,2*expon
1354 goto (10,20,30,30,40) ipot
1355 C----------------------- LJ potential ---------------------------------
1356 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1357 & (sigma0(i),i=1,ntyp)
1359 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1360 write (iout,'(a/)') 'The epsilon array:'
1361 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1362 write (iout,'(/a)') 'One-body parameters:'
1363 write (iout,'(a,4x,a)') 'residue','sigma'
1364 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1367 C----------------------- LJK potential --------------------------------
1368 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1369 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1371 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1372 write (iout,'(a/)') 'The epsilon array:'
1373 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1374 write (iout,'(/a)') 'One-body parameters:'
1375 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1376 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1380 C---------------------- GB or BP potential -----------------------------
1382 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1384 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1385 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1386 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1387 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1388 C now we start reading lipid
1390 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1391 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1393 C print *,"WARNING!!"
1395 C epslip(i,j)=epslip(i,j)+0.05d0
1398 if (lprint) write(iout,*) epslip(1,1),"OK?"
1399 C For the GB potential convert sigma'**2 into chi'
1402 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1406 write (iout,'(/a/)') 'Parameters of the BP potential:'
1407 write (iout,'(a/)') 'The epsilon array:'
1408 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1409 write (iout,'(/a)') 'One-body parameters:'
1410 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1412 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1413 & chip(i),alp(i),i=1,ntyp)
1416 C--------------------- GBV potential -----------------------------------
1417 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1418 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1419 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1421 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1422 write (iout,'(a/)') 'The epsilon array:'
1423 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1424 write (iout,'(/a)') 'One-body parameters:'
1425 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1426 & 's||/s_|_^2',' chip ',' alph '
1427 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1428 & sigii(i),chip(i),alp(i),i=1,ntyp)
1432 C-----------------------------------------------------------------------
1433 C Calculate the "working" parameters of SC interactions.
1437 epslip(i,j)=epslip(j,i)
1442 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1443 sigma(j,i)=sigma(i,j)
1444 rs0(i,j)=dwa16*sigma(i,j)
1448 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1449 & 'Working parameters of the SC interactions:',
1450 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1455 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1464 sigeps=dsign(1.0D0,epsij)
1466 aa_aq(i,j)=epsij*rrij*rrij
1467 bb_aq(i,j)=-sigeps*epsij*rrij
1468 aa_aq(j,i)=aa_aq(i,j)
1469 bb_aq(j,i)=bb_aq(i,j)
1470 epsijlip=epslip(i,j)
1471 sigeps=dsign(1.0D0,epsijlip)
1472 epsijlip=dabs(epsijlip)
1473 aa_lip(i,j)=epsijlip*rrij*rrij
1474 bb_lip(i,j)=-sigeps*epsijlip*rrij
1475 aa_lip(j,i)=aa_lip(i,j)
1476 bb_lip(j,i)=bb_lip(i,j)
1478 sigt1sq=sigma0(i)**2
1479 sigt2sq=sigma0(j)**2
1482 ratsig1=sigt2sq/sigt1sq
1483 ratsig2=1.0D0/ratsig1
1484 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1485 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1486 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1490 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1491 sigmaii(i,j)=rsum_max
1492 sigmaii(j,i)=rsum_max
1494 c sigmaii(i,j)=r0(i,j)
1495 c sigmaii(j,i)=r0(i,j)
1497 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1498 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1499 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1500 augm(i,j)=epsij*r_augm**(2*expon)
1501 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1508 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1509 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1510 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1515 write(iout,*) "tube param"
1516 read(itube,*) epspeptube,sigmapeptube
1517 sigmapeptube=sigmapeptube**6
1518 sigeps=dsign(1.0D0,epspeptube)
1519 epspeptube=dabs(epspeptube)
1520 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1521 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1522 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1524 read(itube,*) epssctube,sigmasctube
1525 sigmasctube=sigmasctube**6
1526 sigeps=dsign(1.0D0,epssctube)
1527 epssctube=dabs(epssctube)
1528 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1529 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1530 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1536 C Define the SC-p interaction constants (hard-coded; old style)
1539 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1541 c aad(i,1)=0.3D0*4.0D0**12
1542 C Following line for constants currently implemented
1543 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1544 aad(i,1)=1.5D0*4.0D0**12
1545 c aad(i,1)=0.17D0*5.6D0**12
1547 C "Soft" SC-p repulsion
1549 C Following line for constants currently implemented
1550 c aad(i,1)=0.3D0*4.0D0**6
1551 C "Hard" SC-p repulsion
1552 bad(i,1)=3.0D0*4.0D0**6
1553 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1562 C 8/9/01 Read the SC-p interaction constants from file
1565 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1568 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1569 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1570 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1571 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1575 write (iout,*) "Parameters of SC-p interactions:"
1577 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1578 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1584 C Define the constants of the disulfide bridge
1588 c Old arbitrary potential - commented out.
1593 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1594 c energy surface of diethyl disulfide.
1595 c A. Liwo and U. Kozlowska, 11/24/03
1611 C if(me.eq.king) then
1612 C write (iout,'(/a)') "Disulfide bridge parameters:"
1613 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1614 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1615 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1616 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1619 if (shield_mode.gt.0) then
1620 C VSolvSphere the volume of solving sphere
1621 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1622 C there will be no distinction between proline peptide group and normal peptide
1623 C group in case of shielding parameters
1624 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
1625 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1626 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1627 write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
1629 C long axis of side chain
1631 long_r_sidechain(i)=vbldsc0(1,i)
1632 short_r_sidechain(i)=sigma0(i)
1637 111 write (iout,*) "Error reading bending energy parameters."
1639 112 write (iout,*) "Error reading rotamer energy parameters."
1641 113 write (iout,*) "Error reading torsional energy parameters."
1643 114 write (iout,*) "Error reading double torsional energy parameters."
1646 & "Error reading cumulant (multibody energy) parameters."
1648 116 write (iout,*) "Error reading electrostatic energy parameters."
1650 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1652 117 write (iout,*) "Error reading side chain interaction parameters."
1654 118 write (iout,*) "Error reading SCp interaction parameters."
1656 119 write (iout,*) "Error reading SCCOR parameters"
1658 121 write (iout,*) "Error in Czybyshev parameters"
1661 call MPI_Finalize(Ierror)
1668 subroutine getenv_loc(var, val)
1669 character(*) var, val
1672 character(2000) line
1675 open (196,file='env',status='old',readonly,shared)
1677 c write(*,*)'looking for ',var
1678 10 read(196,*,err=11,end=11)line
1679 iread=index(line,var)
1680 c write(*,*)iread,' ',var,' ',line
1681 if (iread.eq.0) go to 10
1682 c write(*,*)'---> ',line
1688 iread=iread+ilen(var)+1
1689 read (line(iread:),*,err=12,end=12) val
1690 c write(*,*)'OK: ',var,' = ',val
1696 #elif (defined CRAY)
1697 integer lennam,lenval,ierror
1699 c getenv using a POSIX call, useful on the T3D
1700 c Sept 1996, comment out error check on advice of H. Pritchard
1703 if(lennam.le.0) stop '--error calling getenv--'
1704 call pxfgetenv(var,lennam,val,lenval,ierror)
1705 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1707 call getenv(var,val)
1709 C set the variables used for shielding effect
1710 C if (shield_mode.gt.0) then
1711 C VSolvSphere the volume of solving sphere
1713 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1714 C there will be no distinction between proline peptide group and normal peptide
1715 C group in case of shielding parameters
1716 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1717 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1718 C long axis of side chain
1720 C long_r_sidechain(i)=vbldsc0(1,i)
1721 C short_r_sidechain(i)=sigma0(i)
1723 C lets set the buffor value