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