Merge branch 'master' of mmka.chem.univ.gda.pl:unres
[unres.git] / wham / src-M / parmread.F
1       subroutine parmread(iparm,*)
2 C
3 C Read the parameters of the probability distributions of the virtual-bond
4 C valence angles and the side chains and energy parameters.
5 C
6       implicit real*8 (a-h,o-z)
7       include 'DIMENSIONS'
8       include 'DIMENSIONS.ZSCOPT'
9       include 'DIMENSIONS.FREE'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.CHAIN'
12       include 'COMMON.INTERACT'
13       include 'COMMON.GEO'
14       include 'COMMON.LOCAL'
15       include 'COMMON.TORSION'
16       include 'COMMON.FFIELD'
17       include 'COMMON.NAMES'
18       include 'COMMON.SBRIDGE'
19       include 'COMMON.WEIGHTS'
20       include 'COMMON.ENEPS'
21       include 'COMMON.SCCOR'
22       include 'COMMON.SCROT'
23       include 'COMMON.FREE'
24       include 'COMMON.SHIELD'
25       include 'COMMON.CONTROL'
26       character*1 t1,t2,t3
27       character*1 onelett(4) /"G","A","P","D"/
28       character*1 toronelet(-2:2) /"p","a","G","A","P"/
29       logical lprint
30       dimension blower(3,3,maxlob)
31       character*800 controlcard
32       character*256 bondname_t,thetname_t,rotname_t,torname_t,
33      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
34      & sccorname_t
35       integer ilen
36       external ilen
37       character*16 key
38       integer iparm
39       double precision ip,mp
40       character*6 res1
41 C      write (iout,*) "KURWA"
42 C
43 C Body
44 C
45 C Set LPRINT=.TRUE. for debugging
46       dwa16=2.0d0**(1.0d0/6.0d0)
47       lprint=.false.
48       itypro=20
49 C Assign virtual-bond length
50       vbl=3.8D0
51       vblinv=1.0D0/vbl
52       vblinv2=vblinv*vblinv
53       call card_concat(controlcard,.true.)
54       wname(4)="WCORRH"
55       do i=1,n_ene
56         key = wname(i)(:ilen(wname(i)))
57         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
58       enddo
59
60       write (iout,*) "iparm",iparm," myparm",myparm
61 c If reading not own parameters, skip assignment
62       call reada(controlcard,"D0CM",d0cm,3.78d0)
63       call reada(controlcard,"AKCM",akcm,15.1d0)
64       call reada(controlcard,"AKTH",akth,11.0d0)
65       call reada(controlcard,"AKCT",akct,12.0d0)
66       call reada(controlcard,"V1SS",v1ss,-1.08d0)
67       call reada(controlcard,"V2SS",v2ss,7.61d0)
68       call reada(controlcard,"V3SS",v3ss,13.7d0)
69       call reada(controlcard,"EBR",ebr,-5.50D0)
70       call reada(controlcard,"DTRISS",dtriss,1.0D0)
71       call reada(controlcard,"ATRISS",atriss,0.3D0)
72       call reada(controlcard,"BTRISS",btriss,0.02D0)
73       call reada(controlcard,"CTRISS",ctriss,1.0D0)
74       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
75       write(iout,*) "ATRISS",atriss
76       write(iout,*) "BTRISS",btriss
77       write(iout,*) "CTRISS",ctriss
78       write(iout,*) "DTRISS",dtriss
79
80 C      do i=1,maxres
81 C        dyn_ss_mask(i)=.false.
82 C      enddo
83 C      ebr=-12.0D0
84 c
85 c Old arbitrary potential - commented out.
86 c
87 c      dbr= 4.20D0
88 c      fbr= 3.30D0
89 c
90 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
91 c energy surface of diethyl disulfide.
92 c A. Liwo and U. Kozlowska, 11/24/03
93 c
94       D0CM = 3.78d0
95       AKCM = 15.1d0
96       AKTH = 11.0d0
97       AKCT = 12.0d0
98       V1SS =-1.08d0
99       V2SS = 7.61d0
100       V3SS = 13.7d0
101
102       do i=1,maxres-1
103         do j=i+1,maxres
104           dyn_ssbond_ij(i,j)=1.0d300
105         enddo
106       enddo
107       call reada(controlcard,"HT",Ht,0.0D0)
108 C      if (dyn_ss) then
109 C        ss_depth=ebr/wsc-0.25*eps(1,1)
110 C        write(iout,*) HT,wsc,eps(1,1),'KURWA'
111 C        Ht=Ht/wsc-0.25*eps(1,1)
112        
113 C        akcm=akcm*whpb/wsc
114 C        akth=akth*whpb/wsc
115 C        akct=akct*whpb/wsc
116 C        v1ss=v1ss*whpb/wsc
117 C        v2ss=v2ss*whpb/wsc
118 C        v3ss=v3ss*whpb/wsc
119 C      else
120 C        ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
121 C      endif
122
123       if (iparm.eq.myparm .or. .not.separate_parset) then
124
125 c
126 c Setup weights for UNRES
127 c
128       wsc=ww(1)
129       wscp=ww(2)
130       welec=ww(3)
131       wcorr=ww(4)
132       wcorr5=ww(5)
133       wcorr6=ww(6)
134       wel_loc=ww(7)
135       wturn3=ww(8)
136       wturn4=ww(9)
137       wturn6=ww(10)
138       wang=ww(11)
139       wscloc=ww(12)
140       wtor=ww(13)
141       wtor_d=ww(14)
142       wvdwpp=ww(16)
143       wbond=ww(18)
144       wsccor=ww(19)
145       whpb=ww(15)
146       wstrain=ww(15)
147       wliptran=ww(22)
148       wshield=ww(25)
149       endif
150
151       call card_concat(controlcard,.false.)
152
153 c Return if not own parameters
154
155       if (iparm.ne.myparm .and. separate_parset) return
156
157       call reads(controlcard,"BONDPAR",bondname_t,bondname)
158       open (ibond,file=bondname_t,status='old')
159       rewind(ibond)
160       call reads(controlcard,"THETPAR",thetname_t,thetname)
161       open (ithep,file=thetname_t,status='old')
162       rewind(ithep) 
163       call reads(controlcard,"ROTPAR",rotname_t,rotname)
164       open (irotam,file=rotname_t,status='old')
165       rewind(irotam)
166       call reads(controlcard,"TORPAR",torname_t,torname)
167       open (itorp,file=torname_t,status='old')
168       rewind(itorp)
169       call reads(controlcard,"TORDPAR",tordname_t,tordname)
170       open (itordp,file=tordname_t,status='old')
171       rewind(itordp)
172       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
173       open (isccor,file=sccorname_t,status='old')
174       rewind(isccor)
175       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
176       open (ifourier,file=fouriername_t,status='old')
177       rewind(ifourier)
178       call reads(controlcard,"ELEPAR",elename_t,elename)
179       open (ielep,file=elename_t,status='old')
180       rewind(ielep)
181       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
182       open (isidep,file=sidename_t,status='old')
183       rewind(isidep)
184       call reads(controlcard,"SCPPAR",scpname_t,scpname)
185       open (iscpp,file=scpname_t,status='old')
186       rewind(iscpp)
187       write (iout,*) "Parameter set:",iparm
188       write (iout,*) "Energy-term weights:"
189       do i=1,n_ene
190         write (iout,'(a16,f10.5)') wname(i),ww(i)
191       enddo
192       write (iout,*) "Sidechain potential file        : ",
193      &  sidename_t(:ilen(sidename_t))
194 #ifndef OLDSCP
195       write (iout,*) "SCp potential file              : ",
196      &  scpname_t(:ilen(scpname_t))
197 #endif  
198       write (iout,*) "Electrostatic potential file    : ",
199      &  elename_t(:ilen(elename_t))
200       write (iout,*) "Cumulant coefficient file       : ",
201      &  fouriername_t(:ilen(fouriername_t))
202       write (iout,*) "Torsional parameter file        : ",
203      &  torname_t(:ilen(torname_t))
204       write (iout,*) "Double torsional parameter file : ",
205      &  tordname_t(:ilen(tordname_t))
206       write (iout,*) "Backbone-rotamer parameter file : ",
207      &  sccorname(:ilen(sccorname))
208       write (iout,*) "Bond & inertia constant file    : ",
209      &  bondname_t(:ilen(bondname_t))
210       write (iout,*) "Bending parameter file          : ",
211      &  thetname_t(:ilen(thetname_t))
212       write (iout,*) "Rotamer parameter file          : ",
213      &  rotname_t(:ilen(rotname_t))
214
215 c
216 c Read the virtual-bond parameters, masses, and moments of inertia
217 c and Stokes' radii of the peptide group and side chains
218 c
219 #ifdef CRYST_BOND
220       read (ibond,*) vbldp0,vbldpdum,akp
221       do i=1,ntyp
222         nbondterm(i)=1
223         read (ibond,*) vbldsc0(1,i),aksc(1,i)
224         dsc(i) = vbldsc0(1,i)
225         if (i.eq.10) then
226           dsc_inv(i)=0.0D0
227         else
228           dsc_inv(i)=1.0D0/dsc(i)
229         endif
230       enddo
231 #else
232       read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk
233       do i=1,ntyp
234         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
235      &   j=1,nbondterm(i))
236         dsc(i) = vbldsc0(1,i)
237         if (i.eq.10) then
238           dsc_inv(i)=0.0D0
239         else
240           dsc_inv(i)=1.0D0/dsc(i)
241         endif
242       enddo
243 #endif
244       if (lprint) then
245         write(iout,'(/a/)')"Force constants virtual bonds:"
246         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
247      &   'inertia','Pstok'
248         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
249         do i=1,ntyp
250           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
251      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
252           do j=2,nbondterm(i)
253             write (iout,'(13x,3f10.5)')
254      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
255           enddo
256         enddo
257       endif
258        read(iliptranpar,*) pepliptran
259        do i=1,ntyp
260        read(iliptranpar,*) liptranene(i)
261        enddo
262        close(iliptranpar)
263 #ifdef CRYST_THETA
264 C
265 C Read the parameters of the probability distribution/energy expression 
266 C of the virtual-bond valence angles theta
267 C
268       do i=1,ntyp
269         read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
270      &    (bthet(j,i,1,1),j=1,2)
271         read (ithep,*) (polthet(j,i),j=0,3)
272         read (ithep,*) (gthet(j,i),j=1,3)
273         read (ithep,*) theta0(i),sig0(i),sigc0(i)
274         sigc0(i)=sigc0(i)**2
275       enddo
276       do i=1,ntyp
277       athet(1,i,1,-1)=athet(1,i,1,1)
278       athet(2,i,1,-1)=athet(2,i,1,1)
279       bthet(1,i,1,-1)=-bthet(1,i,1,1)
280       bthet(2,i,1,-1)=-bthet(2,i,1,1)
281       athet(1,i,-1,1)=-athet(1,i,1,1)
282       athet(2,i,-1,1)=-athet(2,i,1,1)
283       bthet(1,i,-1,1)=bthet(1,i,1,1)
284       bthet(2,i,-1,1)=bthet(2,i,1,1)
285       enddo
286       do i=-ntyp,-1
287       a0thet(i)=a0thet(-i)
288       athet(1,i,-1,-1)=athet(1,-i,1,1)
289       athet(2,i,-1,-1)=-athet(2,-i,1,1)
290       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
291       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
292       athet(1,i,-1,1)=athet(1,-i,1,1)
293       athet(2,i,-1,1)=-athet(2,-i,1,1)
294       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
295       bthet(2,i,-1,1)=bthet(2,-i,1,1)
296       athet(1,i,1,-1)=-athet(1,-i,1,1)
297       athet(2,i,1,-1)=athet(2,-i,1,1)
298       bthet(1,i,1,-1)=bthet(1,-i,1,1)
299       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
300       theta0(i)=theta0(-i)
301       sig0(i)=sig0(-i)
302       sigc0(i)=sigc0(-i)
303        do j=0,3
304         polthet(j,i)=polthet(j,-i)
305        enddo
306        do j=1,3
307          gthet(j,i)=gthet(j,-i)
308        enddo
309       enddo
310       close (ithep)
311       if (lprint) then
312 c       write (iout,'(a)') 
313 c    &    'Parameters of the virtual-bond valence angles:'
314 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
315 c    & '    ATHETA0   ','         A1   ','        A2    ',
316 c    & '        B1    ','         B2   '        
317 c       do i=1,ntyp
318 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
319 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
320 c       enddo
321 c       write (iout,'(/a/9x,5a/79(1h-))') 
322 c    & 'Parameters of the expression for sigma(theta_c):',
323 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
324 c    & '     ALPH3    ','    SIGMA0C   '        
325 c       do i=1,ntyp
326 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
327 c    &      (polthet(j,i),j=0,3),sigc0(i) 
328 c       enddo
329 c       write (iout,'(/a/9x,5a/79(1h-))') 
330 c    & 'Parameters of the second gaussian:',
331 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
332 c    & '        G2    ','         G3   '        
333 c       do i=1,ntyp
334 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
335 c    &       sig0(i),(gthet(j,i),j=1,3)
336 c       enddo
337         write (iout,'(a)') 
338      &    'Parameters of the virtual-bond valence angles:'
339         write (iout,'(/a/9x,5a/79(1h-))') 
340      & 'Coefficients of expansion',
341      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
342      & '   b1*10^1    ','    b2*10^1   '        
343         do i=1,ntyp
344           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
345      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
346      &        (10*bthet(j,i,1,1),j=1,2)
347         enddo
348         write (iout,'(/a/9x,5a/79(1h-))') 
349      & 'Parameters of the expression for sigma(theta_c):',
350      & ' alpha0       ','  alph1       ',' alph2        ',
351      & ' alhp3        ','   sigma0c    '        
352         do i=1,ntyp
353           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
354      &      (polthet(j,i),j=0,3),sigc0(i) 
355         enddo
356         write (iout,'(/a/9x,5a/79(1h-))') 
357      & 'Parameters of the second gaussian:',
358      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
359      & '        G2    ','   G3*10^1    '        
360         do i=1,ntyp
361           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
362      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
363         enddo
364       endif
365 #else
366 C
367 C Read the parameters of Utheta determined from ab initio surfaces
368 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
369 C
370       write (iout,*) "tu dochodze"
371       read (ithep,*) nthetyp,ntheterm,ntheterm2,
372      &  ntheterm3,nsingle,ndouble
373       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
374       read (ithep,*) (ithetyp(i),i=1,ntyp1)
375       do i=-ntyp1,-1
376         ithetyp(i)=-ithetyp(-i)
377       enddo
378       write (iout,*) "tu dochodze"
379       do iblock=1,2
380       do i=-maxthetyp,maxthetyp
381         do j=-maxthetyp,maxthetyp
382           do k=-maxthetyp,maxthetyp
383             aa0thet(i,j,k,iblock)=0.0d0
384             do l=1,ntheterm
385               aathet(l,i,j,k,iblock)=0.0d0
386             enddo
387             do l=1,ntheterm2
388               do m=1,nsingle
389                 bbthet(m,l,i,j,k,iblock)=0.0d0
390                 ccthet(m,l,i,j,k,iblock)=0.0d0
391                 ddthet(m,l,i,j,k,iblock)=0.0d0
392                 eethet(m,l,i,j,k,iblock)=0.0d0
393               enddo
394             enddo
395             do l=1,ntheterm3
396               do m=1,ndouble
397                 do mm=1,ndouble
398                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
399                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
400                 enddo
401               enddo
402             enddo
403           enddo
404         enddo
405       enddo
406       enddo
407 C      write (iout,*) "KURWA1"
408       do iblock=1,2
409       do i=0,nthetyp
410         do j=-nthetyp,nthetyp
411           do k=-nthetyp,nthetyp
412             read (ithep,'(6a)') res1
413             write(iout,*) res1,i,j,k
414             read (ithep,*) aa0thet(i,j,k,iblock)
415             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
416             read (ithep,*)
417      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
418      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
419      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
420      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
421      &        ,ll=1,ntheterm2)
422             read (ithep,*)
423      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
424      &      ffthet(lll,llll,ll,i,j,k,iblock),
425      &         ggthet(llll,lll,ll,i,j,k,iblock)
426      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
427      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
428           enddo
429         enddo
430       enddo
431 C       write(iout,*) "KURWA1.1"
432 C
433 C For dummy ends assign glycine-type coefficients of theta-only terms; the
434 C coefficients of theta-and-gamma-dependent terms are zero.
435 C
436       do i=1,nthetyp
437         do j=1,nthetyp
438           do l=1,ntheterm
439             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
440             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
441           enddo
442           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
443           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
444         enddo
445         do l=1,ntheterm
446           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
447         enddo
448         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
449       enddo
450       enddo
451 C       write(iout,*) "KURWA1.5"
452 C Substitution for D aminoacids from symmetry.
453       do iblock=1,2
454       do i=-nthetyp,0
455         do j=-nthetyp,nthetyp
456           do k=-nthetyp,nthetyp
457            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
458            do l=1,ntheterm
459            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
460            enddo
461            do ll=1,ntheterm2
462             do lll=1,nsingle
463             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
464             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
465             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
466             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
467             enddo
468           enddo
469           do ll=1,ntheterm3
470            do lll=2,ndouble
471             do llll=1,lll-1
472             ffthet(llll,lll,ll,i,j,k,iblock)=
473      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
474             ffthet(lll,llll,ll,i,j,k,iblock)=
475      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
476             ggthet(llll,lll,ll,i,j,k,iblock)=
477      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
478             ggthet(lll,llll,ll,i,j,k,iblock)=
479      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
480             enddo !ll
481            enddo  !lll  
482           enddo   !llll
483          enddo    !k
484         enddo     !j
485        enddo      !i
486       enddo       !iblock
487
488 C
489 C Control printout of the coefficients of virtual-bond-angle potentials
490 C
491       if (lprint) then
492         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
493         do i=1,nthetyp+1
494           do j=1,nthetyp+1
495             do k=1,nthetyp+1
496               write (iout,'(//4a)')
497      &         'Type ',onelett(i),onelett(j),onelett(k)
498               write (iout,'(//a,10x,a)') " l","a[l]"
499               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
500               write (iout,'(i2,1pe15.5)')
501      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
502             do l=1,ntheterm2
503               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
504      &          "b",l,"c",l,"d",l,"e",l
505               do m=1,nsingle
506                 write (iout,'(i2,4(1pe15.5))') m,
507      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
508      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
509               enddo
510             enddo
511             do l=1,ntheterm3
512               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
513      &          "f+",l,"f-",l,"g+",l,"g-",l
514               do m=2,ndouble
515                 do n=1,m-1
516                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
517      &              ffthet(n,m,l,i,j,k,iblock),
518      &              ffthet(m,n,l,i,j,k,iblock),
519      &              ggthet(n,m,l,i,j,k,iblock),
520      &              ggthet(m,n,l,i,j,k,iblock)
521                 enddo
522               enddo
523             enddo
524           enddo
525         enddo
526       enddo
527       call flush(iout)
528       endif
529 #endif
530 C      write(iout,*) 'KURWA2'
531 #ifdef CRYST_SC
532 C
533 C Read the parameters of the probability distribution/energy expression
534 C of the side chains.
535 C
536       do i=1,ntyp
537 cc      write (iout,*) "tu dochodze",i
538         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
539         if (i.eq.10) then
540           dsc_inv(i)=0.0D0
541         else
542           dsc_inv(i)=1.0D0/dsc(i)
543         endif
544         if (i.ne.10) then
545         do j=1,nlob(i)
546           do k=1,3
547             do l=1,3
548               blower(l,k,j)=0.0D0
549             enddo
550           enddo
551         enddo  
552         bsc(1,i)=0.0D0
553         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
554         censc(1,1,-i)=censc(1,1,i)
555         censc(2,1,-i)=censc(2,1,i)
556         censc(3,1,-i)=-censc(3,1,i)
557         do j=2,nlob(i)
558           read (irotam,*) bsc(j,i)
559           read (irotam,*) (censc(k,j,i),k=1,3),
560      &                                 ((blower(k,l,j),l=1,k),k=1,3)
561         censc(1,j,-i)=censc(1,j,i)
562         censc(2,j,-i)=censc(2,j,i)
563         censc(3,j,-i)=-censc(3,j,i)
564 C BSC is amplitude of Gaussian
565         enddo
566         do j=1,nlob(i)
567           do k=1,3
568             do l=1,k
569               akl=0.0D0
570               do m=1,3
571                 akl=akl+blower(k,m,j)*blower(l,m,j)
572               enddo
573               gaussc(k,l,j,i)=akl
574               gaussc(l,k,j,i)=akl
575              if (((k.eq.3).and.(l.ne.3))
576      &        .or.((l.eq.3).and.(k.ne.3))) then
577                 gaussc(k,l,j,-i)=-akl
578                 gaussc(l,k,j,-i)=-akl
579               else
580                 gaussc(k,l,j,-i)=akl
581                 gaussc(l,k,j,-i)=akl
582               endif
583             enddo
584           enddo 
585         enddo
586         endif
587       enddo
588       close (irotam)
589       if (lprint) then
590         write (iout,'(/a)') 'Parameters of side-chain local geometry'
591         do i=1,ntyp
592           nlobi=nlob(i)
593           if (nlobi.gt.0) then
594           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
595      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
596 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
597 c          write (iout,'(a,f10.4,4(16x,f10.4))')
598 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
599 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
600            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
601      &                             'log h',(bsc(j,i),j=1,nlobi)
602            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
603      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
604 c          write (iout,'(a)')
605 c         do j=1,nlobi
606 c           ind=0
607 c           do k=1,3
608 c             do l=1,k
609 c              ind=ind+1
610 c              blower(k,l,j)=gaussc(ind,j,i)
611 c             enddo
612 c           enddo
613 c         enddo
614           do k=1,3
615             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
616      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
617           enddo
618           endif
619         enddo
620       endif
621 #else
622 C
623 C Read scrot parameters for potentials determined from all-atom AM1 calculations
624 C added by Urszula Kozlowska 07/11/2007
625 C
626       do i=1,ntyp
627         read (irotam,*)
628        if (i.eq.10) then
629          read (irotam,*)
630        else
631          do j=1,65
632            read(irotam,*) sc_parmin(j,i)
633          enddo
634        endif
635       enddo
636 #endif
637       close(irotam)
638 C      write (iout,*) 'KURWAKURWA'
639 #ifdef CRYST_TOR
640 C
641 C Read torsional parameters in old format
642 C
643       read (itorp,*) ntortyp,nterm_old
644       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
645       read (itorp,*) (itortyp(i),i=1,ntyp)
646       do i=1,ntortyp
647         do j=1,ntortyp
648           read (itorp,'(a)')
649           do k=1,nterm_old
650             read (itorp,*) kk,v1(k,j,i),v2(k,j,i) 
651           enddo
652         enddo
653       enddo
654       close (itorp)
655       if (lprint) then
656         write (iout,'(/a/)') 'Torsional constants:'
657         do i=1,ntortyp
658           do j=1,ntortyp
659             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
660             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
661           enddo
662         enddo
663       endif
664
665
666 #else
667 C
668 C Read torsional parameters
669 C
670       read (itorp,*) ntortyp
671       read (itorp,*) (itortyp(i),i=1,ntyp)
672       write (iout,*) 'ntortyp',ntortyp
673       do iblock=1,2
674       do i=-ntyp,-1
675        itortyp(i)=-itortyp(-i)
676       enddo
677 c      write (iout,*) 'ntortyp',ntortyp
678       do i=0,ntortyp-1
679         do j=-ntortyp+1,ntortyp-1
680           read (itorp,*) nterm(i,j,iblock),
681      &          nlor(i,j,iblock)
682           nterm(-i,-j,iblock)=nterm(i,j,iblock)
683           nlor(-i,-j,iblock)=nlor(i,j,iblock)
684           v0ij=0.0d0
685           si=-1.0d0
686           do k=1,nterm(i,j,iblock)
687             read (itorp,*) kk,v1(k,i,j,iblock),
688      &      v2(k,i,j,iblock)
689             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
690             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
691             v0ij=v0ij+si*v1(k,i,j,iblock)
692             si=-si
693          enddo
694           do k=1,nlor(i,j,iblock)
695             read (itorp,*) kk,vlor1(k,i,j),
696      &        vlor2(k,i,j),vlor3(k,i,j)
697             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
698           enddo
699           v0(i,j,iblock)=v0ij
700           v0(-i,-j,iblock)=v0ij
701         enddo
702       enddo
703       enddo
704       close (itorp)
705       if (lprint) then
706         write (iout,'(/a/)') 'Torsional constants:'
707         do i=1,ntortyp
708           do j=1,ntortyp
709             write (iout,*) 'ityp',i,' jtyp',j
710             write (iout,*) 'Fourier constants'
711             do k=1,nterm(i,j,iblock)
712               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
713      &        v2(k,i,j,iblock)
714             enddo
715             write (iout,*) 'Lorenz constants'
716             do k=1,nlor(i,j,iblock)
717               write (iout,'(3(1pe15.5))')
718      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
719             enddo
720           enddo
721         enddo
722       endif
723 C
724 C 6/23/01 Read parameters for double torsionals
725 C
726       do iblock=1,2
727       do i=0,ntortyp-1
728         do j=-ntortyp+1,ntortyp-1
729           do k=-ntortyp+1,ntortyp-1
730             read (itordp,'(3a1)') t1,t2,t3
731 c              write (iout,*) "OK onelett",
732 c     &         i,j,k,t1,t2,t3
733
734             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
735      &        .or. t3.ne.toronelet(k)) then
736               write (iout,*) "Error in double torsional parameter file",
737      &         i,j,k,t1,t2,t3
738 #ifdef MPI
739               call MPI_Finalize(Ierror)
740 #endif
741                stop "Error in double torsional parameter file"
742             endif
743           read (itordp,*) ntermd_1(i,j,k,iblock),
744      &         ntermd_2(i,j,k,iblock)
745             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
746             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
747             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
748      &         ntermd_1(i,j,k,iblock))
749             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
750      &         ntermd_1(i,j,k,iblock))
751             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
752      &         ntermd_1(i,j,k,iblock))
753             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
754      &         ntermd_1(i,j,k,iblock))
755 C Martix of D parameters for one dimesional foureir series
756             do l=1,ntermd_1(i,j,k,iblock)
757              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
758              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
759              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
760              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
761 c            write(iout,*) "whcodze" ,
762 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
763             enddo
764             read (itordp,*) ((v2c(l,m,i,j,k,iblock),
765      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
766      &         v2s(m,l,i,j,k,iblock),
767      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
768 C Martix of D parameters for two dimesional fourier series
769             do l=1,ntermd_2(i,j,k,iblock)
770              do m=1,l-1
771              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
772              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
773              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
774              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
775              enddo!m
776             enddo!l
777           enddo!k
778         enddo!j
779       enddo!i
780       enddo!iblock
781       if (lprint) then
782       write (iout,*)
783       write (iout,*) 'Constants for double torsionals'
784       do iblock=1,2
785       do i=0,ntortyp-1
786         do j=-ntortyp+1,ntortyp-1
787           do k=-ntortyp+1,ntortyp-1
788             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
789      &        ' nsingle',ntermd_1(i,j,k,iblock),
790      &        ' ndouble',ntermd_2(i,j,k,iblock)
791             write (iout,*)
792             write (iout,*) 'Single angles:'
793             do l=1,ntermd_1(i,j,k,iblock)
794               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
795      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
796      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
797      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
798             enddo
799             write (iout,*)
800             write (iout,*) 'Pairs of angles:'
801             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
802             do l=1,ntermd_2(i,j,k,iblock)
803               write (iout,'(i5,20f10.5)')
804      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
805             enddo
806             write (iout,*)
807            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
808             do l=1,ntermd_2(i,j,k,iblock)
809               write (iout,'(i5,20f10.5)')
810      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
811      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
812             enddo
813             write (iout,*)
814           enddo
815         enddo
816       enddo
817       enddo
818       endif
819 #endif
820 C Read of Side-chain backbone correlation parameters
821 C Modified 11 May 2012 by Adasko
822 CCC
823 C
824       read (isccor,*) nsccortyp
825       read (isccor,*) (isccortyp(i),i=1,ntyp)
826       do i=-ntyp,-1
827         isccortyp(i)=-isccortyp(-i)
828       enddo
829       iscprol=isccortyp(20)
830 c      write (iout,*) 'ntortyp',ntortyp
831       maxinter=3
832 cc maxinter is maximum interaction sites
833       do l=1,maxinter
834       do i=1,nsccortyp
835         do j=1,nsccortyp
836           read (isccor,*)
837      &nterm_sccor(i,j),nlor_sccor(i,j)
838           write (iout,*) nterm_sccor(i,j)
839           v0ijsccor=0.0d0
840           v0ijsccor1=0.0d0
841           v0ijsccor2=0.0d0
842           v0ijsccor3=0.0d0
843           si=-1.0d0
844           nterm_sccor(-i,j)=nterm_sccor(i,j)
845           nterm_sccor(-i,-j)=nterm_sccor(i,j)
846           nterm_sccor(i,-j)=nterm_sccor(i,j)
847           write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
848      &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
849           do k=1,nterm_sccor(i,j)
850             read (isccor,*) kk,v1sccor(k,l,i,j)
851      &    ,v2sccor(k,l,i,j)
852             if (j.eq.iscprol) then
853              if (i.eq.isccortyp(10)) then
854              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
855              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
856              else
857              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
858      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
859              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
860      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
861              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
862              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
863              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
864              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
865              endif
866             else
867              if (i.eq.isccortyp(10)) then
868              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
869              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
870              else
871                if (j.eq.isccortyp(10)) then
872              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
873              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
874                else
875              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
876              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
877              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
878              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
879              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
880              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
881                 endif
882                endif
883             endif
884             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
885             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
886             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
887             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
888             si=-si
889            enddo
890           do k=1,nlor_sccor(i,j)
891             read (isccor,*) kk,vlor1sccor(k,i,j),
892      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
893             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
894      &(1+vlor3sccor(k,i,j)**2)
895           enddo
896           v0sccor(l,i,j)=v0ijsccor
897           v0sccor(l,-i,j)=v0ijsccor1
898           v0sccor(l,i,-j)=v0ijsccor2
899           v0sccor(l,-i,-j)=v0ijsccor3
900           enddo
901         enddo
902       enddo
903       close (isccor)
904       if (lprint) then
905         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
906         do i=1,nsccortyp
907           do j=1,nsccortyp
908             write (iout,*) 'ityp',i,' jtyp',j
909             write (iout,*) 'Fourier constants'
910             do k=1,nterm_sccor(i,j)
911               write (iout,'(2(1pe15.5))')
912      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
913             enddo
914             write (iout,*) 'Lorenz constants'
915             do k=1,nlor_sccor(i,j)
916               write (iout,'(3(1pe15.5))')
917      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
918             enddo
919           enddo
920         enddo
921       endif
922 C
923 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
924 C         interaction energy of the Gly, Ala, and Pro prototypes.
925 C
926       read (ifourier,*) nloctyp
927       do i=0,nloctyp-1
928         read (ifourier,*)
929         read (ifourier,*) (b(ii,i),ii=1,13)
930         if (lprint) then
931         write (iout,*) 'Type',i
932         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
933         endif
934         B1(1,i)  = b(3,i)
935         B1(2,i)  = b(5,i)
936         B1(1,-i) = b(3,i)
937         B1(2,-i) = -b(5,i)
938 c        b1(1,i)=0.0d0
939 c        b1(2,i)=0.0d0
940         B1tilde(1,i) = b(3,i)
941         B1tilde(2,i) =-b(5,i)
942         B1tilde(1,-i) =-b(3,i)
943         B1tilde(2,-i) =b(5,i)
944 c        b1tilde(1,i)=0.0d0
945 c        b1tilde(2,i)=0.0d0
946         B2(1,i)  = b(2,i)
947         B2(2,i)  = b(4,i)
948         B2(1,-i)  =b(2,i)
949         B2(2,-i)  =-b(4,i)
950
951 c        b2(1,i)=0.0d0
952 c        b2(2,i)=0.0d0
953         CC(1,1,i)= b(7,i)
954         CC(2,2,i)=-b(7,i)
955         CC(2,1,i)= b(9,i)
956         CC(1,2,i)= b(9,i)
957         CC(1,1,-i)= b(7,i)
958         CC(2,2,-i)=-b(7,i)
959         CC(2,1,-i)=-b(9,i)
960         CC(1,2,-i)=-b(9,i)
961 c        CC(1,1,i)=0.0d0
962 c        CC(2,2,i)=0.0d0
963 c        CC(2,1,i)=0.0d0
964 c        CC(1,2,i)=0.0d0
965         Ctilde(1,1,i)=b(7,i)
966         Ctilde(1,2,i)=b(9,i)
967         Ctilde(2,1,i)=-b(9,i)
968         Ctilde(2,2,i)=b(7,i)
969         Ctilde(1,1,-i)=b(7,i)
970         Ctilde(1,2,-i)=-b(9,i)
971         Ctilde(2,1,-i)=b(9,i)
972         Ctilde(2,2,-i)=b(7,i)
973
974 c        Ctilde(1,1,i)=0.0d0
975 c        Ctilde(1,2,i)=0.0d0
976 c        Ctilde(2,1,i)=0.0d0
977 c        Ctilde(2,2,i)=0.0d0
978         DD(1,1,i)= b(6,i)
979         DD(2,2,i)=-b(6,i)
980         DD(2,1,i)= b(8,i)
981         DD(1,2,i)= b(8,i)
982         DD(1,1,-i)= b(6,i)
983         DD(2,2,-i)=-b(6,i)
984         DD(2,1,-i)=-b(8,i)
985         DD(1,2,-i)=-b(8,i)
986 c        DD(1,1,i)=0.0d0
987 c        DD(2,2,i)=0.0d0
988 c        DD(2,1,i)=0.0d0
989 c        DD(1,2,i)=0.0d0
990         Dtilde(1,1,i)=b(6,i)
991         Dtilde(1,2,i)=b(8,i)
992         Dtilde(2,1,i)=-b(8,i)
993         Dtilde(2,2,i)=b(6,i)
994         Dtilde(1,1,-i)=b(6,i)
995         Dtilde(1,2,-i)=-b(8,i)
996         Dtilde(2,1,-i)=b(8,i)
997         Dtilde(2,2,-i)=b(6,i)
998
999 c        Dtilde(1,1,i)=0.0d0
1000 c        Dtilde(1,2,i)=0.0d0
1001 c        Dtilde(2,1,i)=0.0d0
1002 c        Dtilde(2,2,i)=0.0d0
1003         EE(1,1,i)= b(10,i)+b(11,i)
1004         EE(2,2,i)=-b(10,i)+b(11,i)
1005         EE(2,1,i)= b(12,i)-b(13,i)
1006         EE(1,2,i)= b(12,i)+b(13,i)
1007         EE(1,1,-i)= b(10,i)+b(11,i)
1008         EE(2,2,-i)=-b(10,i)+b(11,i)
1009         EE(2,1,-i)=-b(12,i)+b(13,i)
1010         EE(1,2,-i)=-b(12,i)-b(13,i)
1011
1012 c        ee(1,1,i)=1.0d0
1013 c        ee(2,2,i)=1.0d0
1014 c        ee(2,1,i)=0.0d0
1015 c        ee(1,2,i)=0.0d0
1016 c        ee(2,1,i)=ee(1,2,i)
1017
1018       enddo
1019       if (lprint) then
1020       do i=1,nloctyp
1021         write (iout,*) 'Type',i
1022         write (iout,*) 'B1'
1023 c        write (iout,'(f10.5)') B1(:,i)
1024         write(iout,*) B1(1,i),B1(2,i)
1025         write (iout,*) 'B2'
1026 c        write (iout,'(f10.5)') B2(:,i)
1027         write(iout,*) B2(1,i),B2(2,i)
1028         write (iout,*) 'CC'
1029         do j=1,2
1030           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1031         enddo
1032         write(iout,*) 'DD'
1033         do j=1,2
1034           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1035         enddo
1036         write(iout,*) 'EE'
1037         do j=1,2
1038           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1039         enddo
1040       enddo
1041       endif
1042
1043 C Read electrostatic-interaction parameters
1044 C
1045       if (lprint) then
1046         write (iout,'(/a)') 'Electrostatic interaction constants:'
1047         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
1048      &            'IT','JT','APP','BPP','AEL6','AEL3'
1049       endif
1050       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1051       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1052       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1053       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1054       close (ielep)
1055       do i=1,2
1056         do j=1,2
1057         rri=rpp(i,j)**6
1058         app (i,j)=epp(i,j)*rri*rri 
1059         bpp (i,j)=-2.0D0*epp(i,j)*rri
1060         ael6(i,j)=elpp6(i,j)*4.2D0**6
1061         ael3(i,j)=elpp3(i,j)*4.2D0**3
1062         lprint=.true.
1063         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1064      &                    ael6(i,j),ael3(i,j)
1065          lprint=.false.
1066         enddo
1067       enddo
1068 C
1069 C Read side-chain interaction parameters.
1070 C
1071       read (isidep,*) ipot,expon
1072       if (ipot.lt.1 .or. ipot.gt.5) then
1073         write (iout,'(2a)') 'Error while reading SC interaction',
1074      &               'potential file - unknown potential type.'
1075         stop
1076       endif
1077       expon2=expon/2
1078       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1079      & ', exponents are ',expon,2*expon 
1080       goto (10,20,30,30,40) ipot
1081 C----------------------- LJ potential ---------------------------------
1082    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1083       if (lprint) then
1084         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1085         write (iout,'(a/)') 'The epsilon array:'
1086         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1087         write (iout,'(/a)') 'One-body parameters:'
1088         write (iout,'(a,4x,a)') 'residue','sigma'
1089         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1090       endif
1091       goto 50
1092 C----------------------- LJK potential --------------------------------
1093    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1094      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1095       if (lprint) then
1096         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1097         write (iout,'(a/)') 'The epsilon array:'
1098         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1099         write (iout,'(/a)') 'One-body parameters:'
1100         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1101         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1102      &        i=1,ntyp)
1103       endif
1104       goto 50
1105 C---------------------- GB or BP potential -----------------------------
1106    30 do i=1,ntyp
1107        read (isidep,*)(eps(i,j),j=i,ntyp)
1108       enddo
1109       read (isidep,*)(sigma0(i),i=1,ntyp)
1110       read (isidep,*)(sigii(i),i=1,ntyp)
1111       read (isidep,*)(chip(i),i=1,ntyp)
1112       read (isidep,*)(alp(i),i=1,ntyp)
1113       do i=1,ntyp
1114        read (isidep,*)(epslip(i,j),j=i,ntyp)
1115 C       write(iout,*) "WARNING!!",i,ntyp
1116        write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1117 C       do j=1,ntyp
1118 C       epslip(i,j)=epslip(i,j)+0.05d0
1119 C       enddo
1120       enddo
1121 C For the GB potential convert sigma'**2 into chi'
1122       if (ipot.eq.4) then
1123         do i=1,ntyp
1124           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1125         enddo
1126       endif
1127       if (lprint) then
1128         write (iout,'(/a/)') 'Parameters of the BP potential:'
1129         write (iout,'(a/)') 'The epsilon array:'
1130         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1131         write (iout,'(/a)') 'One-body parameters:'
1132         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1133      &       '    chip  ','    alph  '
1134         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1135      &                     chip(i),alp(i),i=1,ntyp)
1136       endif
1137       goto 50
1138 C--------------------- GBV potential -----------------------------------
1139    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1140      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1141      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1142       if (lprint) then
1143         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1144         write (iout,'(a/)') 'The epsilon array:'
1145         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1146         write (iout,'(/a)') 'One-body parameters:'
1147         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1148      &      's||/s_|_^2','    chip  ','    alph  '
1149         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1150      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1151       endif
1152    50 continue
1153       close (isidep)
1154 C-----------------------------------------------------------------------
1155 C Calculate the "working" parameters of SC interactions.
1156       do i=2,ntyp
1157         do j=1,i-1
1158           eps(i,j)=eps(j,i)
1159           epslip(i,j)=epslip(j,i)
1160         enddo
1161       enddo
1162       do i=1,ntyp
1163         do j=i,ntyp
1164           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1165           sigma(j,i)=sigma(i,j)
1166           rs0(i,j)=dwa16*sigma(i,j)
1167           rs0(j,i)=rs0(i,j)
1168         enddo
1169       enddo
1170       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
1171      & 'Working parameters of the SC interactions:',
1172      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1173      & '  chi1   ','   chi2   ' 
1174       do i=1,ntyp
1175         do j=i,ntyp
1176           epsij=eps(i,j)
1177           epsijlip=epslip(i,j)
1178           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1179             rrij=sigma(i,j)
1180           else
1181             rrij=rr0(i)+rr0(j)
1182           endif
1183           r0(i,j)=rrij
1184           r0(j,i)=rrij
1185           rrij=rrij**expon
1186           epsij=eps(i,j)
1187           sigeps=dsign(1.0D0,epsij)
1188           epsij=dabs(epsij)
1189           aa_aq(i,j)=epsij*rrij*rrij
1190           bb_aq(i,j)=-sigeps*epsij*rrij
1191           aa_aq(j,i)=aa_aq(i,j)
1192           bb_aq(j,i)=bb_aq(i,j)
1193           sigeps=dsign(1.0D0,epsijlip)
1194           epsijlip=dabs(epsijlip)
1195           aa_lip(i,j)=epsijlip*rrij*rrij
1196           bb_lip(i,j)=-sigeps*epsijlip*rrij
1197           aa_lip(j,i)=aa_lip(i,j)
1198           bb_lip(j,i)=bb_lip(i,j)
1199           if (ipot.gt.2) then
1200             sigt1sq=sigma0(i)**2
1201             sigt2sq=sigma0(j)**2
1202             sigii1=sigii(i)
1203             sigii2=sigii(j)
1204             ratsig1=sigt2sq/sigt1sq
1205             ratsig2=1.0D0/ratsig1
1206             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1207             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1208             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1209           else
1210             rsum_max=sigma(i,j)
1211           endif
1212 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1213             sigmaii(i,j)=rsum_max
1214             sigmaii(j,i)=rsum_max 
1215 c         else
1216 c           sigmaii(i,j)=r0(i,j)
1217 c           sigmaii(j,i)=r0(i,j)
1218 c         endif
1219 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1220           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1221             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1222             augm(i,j)=epsij*r_augm**(2*expon)
1223 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1224             augm(j,i)=augm(i,j)
1225           else
1226             augm(i,j)=0.0D0
1227             augm(j,i)=0.0D0
1228           endif
1229           if (lprint) then
1230             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
1231      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1232      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1233           endif
1234         enddo
1235       enddo
1236 C
1237 C Define the SC-p interaction constants
1238 C
1239 #ifdef OLDSCP
1240       do i=1,20
1241 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
1242 C helix formation)
1243 c       aad(i,1)=0.3D0*4.0D0**12
1244 C Following line for constants currently implemented
1245 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1246         aad(i,1)=1.5D0*4.0D0**12
1247 c       aad(i,1)=0.17D0*5.6D0**12
1248         aad(i,2)=aad(i,1)
1249 C "Soft" SC-p repulsion
1250         bad(i,1)=0.0D0
1251 C Following line for constants currently implemented
1252 c       aad(i,1)=0.3D0*4.0D0**6
1253 C "Hard" SC-p repulsion
1254         bad(i,1)=3.0D0*4.0D0**6
1255 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1256         bad(i,2)=bad(i,1)
1257 c       aad(i,1)=0.0D0
1258 c       aad(i,2)=0.0D0
1259 c       bad(i,1)=1228.8D0
1260 c       bad(i,2)=1228.8D0
1261       enddo
1262 #else
1263 C
1264 C 8/9/01 Read the SC-p interaction constants from file
1265 C
1266       do i=1,ntyp
1267         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1268       enddo
1269       do i=1,ntyp
1270         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1271         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1272         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1273         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1274       enddo
1275
1276       if (lprint) then
1277         write (iout,*) "Parameters of SC-p interactions:"
1278         do i=1,20
1279           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1280      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1281         enddo
1282       endif
1283 #endif
1284 C
1285 C Define the constants of the disulfide bridge
1286 C
1287 C      ebr=-12.0D0
1288 c
1289 c Old arbitrary potential - commented out.
1290 c
1291 c      dbr= 4.20D0
1292 c      fbr= 3.30D0
1293 c
1294 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1295 c energy surface of diethyl disulfide.
1296 c A. Liwo and U. Kozlowska, 11/24/03
1297 c
1298 C      D0CM = 3.78d0
1299 C      AKCM = 15.1d0
1300 C      AKTH = 11.0d0
1301 C      AKCT = 12.0d0
1302 C      V1SS =-1.08d0
1303 C      V2SS = 7.61d0
1304 C      V3SS = 13.7d0
1305       write (iout,*) dyn_ss,'dyndyn'
1306       if (dyn_ss) then
1307         ss_depth=ebr/wsc-0.25*eps(1,1)
1308 C        write(iout,*) akcm,whpb,wsc,'KURWA'
1309         Ht=Ht/wsc-0.25*eps(1,1)
1310
1311         akcm=akcm*whpb/wsc
1312         akth=akth*whpb/wsc
1313         akct=akct*whpb/wsc
1314         v1ss=v1ss*whpb/wsc
1315         v2ss=v2ss*whpb/wsc
1316         v3ss=v3ss*whpb/wsc
1317       else
1318         ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1319       endif
1320
1321 C      if (lprint) then
1322       write (iout,'(/a)') "Disulfide bridge parameters:"
1323       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1324       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1325       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1326       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1327      & ' v3ss:',v3ss
1328 C      endif
1329       if (shield_mode.gt.0) then
1330       pi=3.141592d0
1331 C VSolvSphere the volume of solving sphere
1332 C      print *,pi,"pi"
1333 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
1334 C there will be no distinction between proline peptide group and normal peptide
1335 C group in case of shielding parameters
1336       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1337       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1338       write (iout,*) VSolvSphere,VSolvSphere_div
1339 C long axis of side chain 
1340       do i=1,ntyp
1341       long_r_sidechain(i)=vbldsc0(1,i)
1342       short_r_sidechain(i)=sigma0(i)
1343       enddo
1344       buff_shield=1.0d0
1345       endif 
1346       return
1347       end