8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
871 read(iion,*) msc(i,5),restok(i,5)
872 print *,msc(i,5),restok(i,5)
876 read (iion,*) ncatprotparm
877 allocate(catprm(ncatprotparm,4))
879 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
882 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
883 !----------------------------------------------------
884 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
885 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
886 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
887 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
888 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
889 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
899 allocate(liptranene(ntyp))
900 !C reading lipid parameters
901 write (iout,*) "iliptranpar",iliptranpar
903 read(iliptranpar,*) pepliptran
906 read(iliptranpar,*) liptranene(i)
907 print *,liptranene(i)
913 ! Read the parameters of the probability distribution/energy expression
914 ! of the virtual-bond valence angles theta
917 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
918 (bthet(j,i,1,1),j=1,2)
919 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
920 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
921 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
925 athet(1,i,1,-1)=athet(1,i,1,1)
926 athet(2,i,1,-1)=athet(2,i,1,1)
927 bthet(1,i,1,-1)=-bthet(1,i,1,1)
928 bthet(2,i,1,-1)=-bthet(2,i,1,1)
929 athet(1,i,-1,1)=-athet(1,i,1,1)
930 athet(2,i,-1,1)=-athet(2,i,1,1)
931 bthet(1,i,-1,1)=bthet(1,i,1,1)
932 bthet(2,i,-1,1)=bthet(2,i,1,1)
936 athet(1,i,-1,-1)=athet(1,-i,1,1)
937 athet(2,i,-1,-1)=-athet(2,-i,1,1)
938 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
939 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
940 athet(1,i,-1,1)=athet(1,-i,1,1)
941 athet(2,i,-1,1)=-athet(2,-i,1,1)
942 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
943 bthet(2,i,-1,1)=bthet(2,-i,1,1)
944 athet(1,i,1,-1)=-athet(1,-i,1,1)
945 athet(2,i,1,-1)=athet(2,-i,1,1)
946 bthet(1,i,1,-1)=bthet(1,-i,1,1)
947 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
952 polthet(j,i)=polthet(j,-i)
955 gthet(j,i)=gthet(j,-i)
963 'Parameters of the virtual-bond valence angles:'
964 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
965 ' ATHETA0 ',' A1 ',' A2 ',&
968 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
969 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
971 write (iout,'(/a/9x,5a/79(1h-))') &
972 'Parameters of the expression for sigma(theta_c):',&
973 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
974 ' ALPH3 ',' SIGMA0C '
976 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
977 (polthet(j,i),j=0,3),sigc0(i)
979 write (iout,'(/a/9x,5a/79(1h-))') &
980 'Parameters of the second gaussian:',&
981 ' THETA0 ',' SIGMA0 ',' G1 ',&
984 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
985 sig0(i),(gthet(j,i),j=1,3)
989 'Parameters of the virtual-bond valence angles:'
990 write (iout,'(/a/9x,5a/79(1h-))') &
991 'Coefficients of expansion',&
992 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
993 ' b1*10^1 ',' b2*10^1 '
995 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
996 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
997 (10*bthet(j,i,1,1),j=1,2)
999 write (iout,'(/a/9x,5a/79(1h-))') &
1000 'Parameters of the expression for sigma(theta_c):',&
1001 ' alpha0 ',' alph1 ',' alph2 ',&
1002 ' alhp3 ',' sigma0c '
1004 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1005 (polthet(j,i),j=0,3),sigc0(i)
1007 write (iout,'(/a/9x,5a/79(1h-))') &
1008 'Parameters of the second gaussian:',&
1009 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1012 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1013 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1019 ! Read the parameters of Utheta determined from ab initio surfaces
1020 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1022 IF (tor_mode.eq.0) THEN
1023 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1024 ntheterm3,nsingle,ndouble
1025 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1027 !----------------------------------------------------
1028 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1029 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1030 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1031 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1032 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1033 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1034 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1035 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1036 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1037 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1038 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1039 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1040 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1041 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1042 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1043 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1044 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1045 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1046 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1047 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1048 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1049 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1050 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1051 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1053 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1055 ithetyp(i)=-ithetyp(-i)
1058 aa0thet(:,:,:,:)=0.0d0
1059 aathet(:,:,:,:,:)=0.0d0
1060 bbthet(:,:,:,:,:,:)=0.0d0
1061 ccthet(:,:,:,:,:,:)=0.0d0
1062 ddthet(:,:,:,:,:,:)=0.0d0
1063 eethet(:,:,:,:,:,:)=0.0d0
1064 ffthet(:,:,:,:,:,:,:)=0.0d0
1065 ggthet(:,:,:,:,:,:,:)=0.0d0
1067 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1069 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1070 ! VAR:1=non-glicyne non-proline 2=proline
1071 ! VAR:negative values for D-aminoacid
1073 do j=-nthetyp,nthetyp
1074 do k=-nthetyp,nthetyp
1075 read (ithep,'(6a)',end=111,err=111) res1
1076 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1077 ! VAR: aa0thet is variable describing the average value of Foureir
1078 ! VAR: expansion series
1079 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1080 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1081 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1082 read (ithep,*,end=111,err=111) &
1083 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1084 read (ithep,*,end=111,err=111) &
1085 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1086 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1087 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1088 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1090 read (ithep,*,end=111,err=111) &
1091 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1092 ffthet(lll,llll,ll,i,j,k,iblock),&
1093 ggthet(llll,lll,ll,i,j,k,iblock),&
1094 ggthet(lll,llll,ll,i,j,k,iblock),&
1095 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1100 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1101 ! coefficients of theta-and-gamma-dependent terms are zero.
1102 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1103 ! RECOMENTDED AFTER VERSION 3.3)
1107 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1108 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1110 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1111 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1114 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1116 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1119 ! AND COMMENT THE LOOPS BELOW
1123 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1124 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1126 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1127 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1130 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1132 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1137 ! Substitution for D aminoacids from symmetry.
1140 do j=-nthetyp,nthetyp
1141 do k=-nthetyp,nthetyp
1142 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1144 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1148 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1149 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1150 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1151 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1157 ffthet(llll,lll,ll,i,j,k,iblock)= &
1158 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1159 ffthet(lll,llll,ll,i,j,k,iblock)= &
1160 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1161 ggthet(llll,lll,ll,i,j,k,iblock)= &
1162 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1163 ggthet(lll,llll,ll,i,j,k,iblock)= &
1164 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1173 ! Control printout of the coefficients of virtual-bond-angle potentials
1176 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1181 write (iout,'(//4a)') &
1182 'Type ',onelett(i),onelett(j),onelett(k)
1183 write (iout,'(//a,10x,a)') " l","a[l]"
1184 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1185 write (iout,'(i2,1pe15.5)') &
1186 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1188 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1189 "b",l,"c",l,"d",l,"e",l
1191 write (iout,'(i2,4(1pe15.5))') m,&
1192 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1193 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1197 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1198 "f+",l,"f-",l,"g+",l,"g-",l
1201 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1202 ffthet(n,m,l,i,j,k,iblock),&
1203 ffthet(m,n,l,i,j,k,iblock),&
1204 ggthet(n,m,l,i,j,k,iblock),&
1205 ggthet(m,n,l,i,j,k,iblock)
1216 !C here will be the apropriate recalibrating for D-aminoacid
1217 read (ithep,*,end=121,err=121) nthetyp
1218 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1219 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1220 do i=-nthetyp+1,nthetyp-1
1221 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1222 do j=0,nbend_kcc_Tb(i)
1223 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1227 write (iout,'(a)') &
1228 "Parameters of the valence-only potentials"
1229 do i=-nthetyp+1,nthetyp-1
1230 write (iout,'(2a)') "Type ",toronelet(i)
1231 do k=0,nbend_kcc_Tb(i)
1232 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1238 write (2,*) "Start reading THETA_PDB",ithep_pdb
1240 ! write (2,*) 'i=',i
1241 read (ithep_pdb,*,err=111,end=111) &
1242 a0thet(i),(athet(j,i,1,1),j=1,2),&
1243 (bthet(j,i,1,1),j=1,2)
1244 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1245 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1246 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1247 sigc0(i)=sigc0(i)**2
1250 athet(1,i,1,-1)=athet(1,i,1,1)
1251 athet(2,i,1,-1)=athet(2,i,1,1)
1252 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1253 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1254 athet(1,i,-1,1)=-athet(1,i,1,1)
1255 athet(2,i,-1,1)=-athet(2,i,1,1)
1256 bthet(1,i,-1,1)=bthet(1,i,1,1)
1257 bthet(2,i,-1,1)=bthet(2,i,1,1)
1260 a0thet(i)=a0thet(-i)
1261 athet(1,i,-1,-1)=athet(1,-i,1,1)
1262 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1263 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1264 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1265 athet(1,i,-1,1)=athet(1,-i,1,1)
1266 athet(2,i,-1,1)=-athet(2,-i,1,1)
1267 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1268 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1269 athet(1,i,1,-1)=-athet(1,-i,1,1)
1270 athet(2,i,1,-1)=athet(2,-i,1,1)
1271 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1272 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1273 theta0(i)=theta0(-i)
1277 polthet(j,i)=polthet(j,-i)
1280 gthet(j,i)=gthet(j,-i)
1283 write (2,*) "End reading THETA_PDB"
1287 !--------------- Reading theta parameters for nucleic acid-------
1288 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1289 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1290 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1291 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1292 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1293 nthetyp_nucl+1,nthetyp_nucl+1))
1294 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1295 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1296 nthetyp_nucl+1,nthetyp_nucl+1))
1297 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1298 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1299 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1300 nthetyp_nucl+1,nthetyp_nucl+1))
1301 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1302 nthetyp_nucl+1,nthetyp_nucl+1))
1303 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1304 nthetyp_nucl+1,nthetyp_nucl+1))
1305 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1306 nthetyp_nucl+1,nthetyp_nucl+1))
1307 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1308 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1309 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1310 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1311 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1312 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1314 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1315 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1317 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1319 aa0thet_nucl(:,:,:)=0.0d0
1320 aathet_nucl(:,:,:,:)=0.0d0
1321 bbthet_nucl(:,:,:,:,:)=0.0d0
1322 ccthet_nucl(:,:,:,:,:)=0.0d0
1323 ddthet_nucl(:,:,:,:,:)=0.0d0
1324 eethet_nucl(:,:,:,:,:)=0.0d0
1325 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1326 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1331 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1332 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1333 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1334 read (ithep_nucl,*,end=111,err=111) &
1335 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1336 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1337 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1338 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1339 read (ithep_nucl,*,end=111,err=111) &
1340 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1341 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1342 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1347 !-------------------------------------------
1348 allocate(nlob(ntyp1)) !(ntyp1)
1349 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1350 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1351 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1358 gaussc(:,:,:,:)=0.0D0
1362 ! Read the parameters of the probability distribution/energy expression
1363 ! of the side chains.
1366 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1370 dsc_inv(i)=1.0D0/dsc(i)
1381 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1382 ((blower(k,l,1),l=1,k),k=1,3)
1383 censc(1,1,-i)=censc(1,1,i)
1384 censc(2,1,-i)=censc(2,1,i)
1385 censc(3,1,-i)=-censc(3,1,i)
1387 read (irotam,*,end=112,err=112) bsc(j,i)
1388 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1389 ((blower(k,l,j),l=1,k),k=1,3)
1390 censc(1,j,-i)=censc(1,j,i)
1391 censc(2,j,-i)=censc(2,j,i)
1392 censc(3,j,-i)=-censc(3,j,i)
1393 ! BSC is amplitude of Gaussian
1400 akl=akl+blower(k,m,j)*blower(l,m,j)
1404 if (((k.eq.3).and.(l.ne.3)) &
1405 .or.((l.eq.3).and.(k.ne.3))) then
1406 gaussc(k,l,j,-i)=-akl
1407 gaussc(l,k,j,-i)=-akl
1409 gaussc(k,l,j,-i)=akl
1410 gaussc(l,k,j,-i)=akl
1419 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1422 if (nlobi.gt.0) then
1424 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1425 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1426 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1427 'log h',(bsc(j,i),j=1,nlobi)
1428 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1429 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1431 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1432 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1435 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1436 write (iout,'(a,f10.4,4(16x,f10.4))') &
1437 'Center ',(bsc(j,i),j=1,nlobi)
1438 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1447 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1448 ! added by Urszula Kozlowska 07/11/2007
1450 !el Maximum number of SC local term fitting function coefficiants
1451 !el integer,parameter :: maxsccoef=65
1453 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1456 read (irotam,*,end=112,err=112)
1458 read (irotam,*,end=112,err=112)
1461 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1465 !---------reading nucleic acid parameters for rotamers-------------------
1466 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1467 do i=1,ntyp_molec(2)
1468 read (irotam_nucl,*,end=112,err=112)
1470 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1476 write (iout,*) "Base rotamer parameters"
1477 do i=1,ntyp_molec(2)
1478 write (iout,'(a)') restyp(i,2)
1479 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1484 ! Read the parameters of the probability distribution/energy expression
1485 ! of the side chains.
1487 write (2,*) "Start reading ROTAM_PDB"
1489 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1493 dsc_inv(i)=1.0D0/dsc(i)
1504 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1505 ((blower(k,l,1),l=1,k),k=1,3)
1507 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1508 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1509 ((blower(k,l,j),l=1,k),k=1,3)
1516 akl=akl+blower(k,m,j)*blower(l,m,j)
1526 write (2,*) "End reading ROTAM_PDB"
1532 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1533 !C interaction energy of the Gly, Ala, and Pro prototypes.
1535 read (ifourier,*) nloctyp
1536 SPLIT_FOURIERTOR = nloctyp.lt.0
1537 nloctyp = iabs(nloctyp)
1538 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1539 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1540 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1541 !C allocate(ctilde(2,2,nres))
1542 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1543 !C allocate(gtb1(2,nres))
1544 !C allocate(gtb2(2,nres))
1545 !C allocate(cc(2,2,nres))
1546 !C allocate(dd(2,2,nres))
1547 !C allocate(ee(2,2,nres))
1548 !C allocate(gtcc(2,2,nres))
1549 !C allocate(gtdd(2,2,nres))
1550 !C allocate(gtee(2,2,nres))
1553 allocate(itype2loc(-ntyp1:ntyp1))
1554 allocate(iloctyp(-nloctyp:nloctyp))
1555 allocate(bnew1(3,2,-nloctyp:nloctyp))
1556 allocate(bnew2(3,2,-nloctyp:nloctyp))
1557 allocate(ccnew(3,2,-nloctyp:nloctyp))
1558 allocate(ddnew(3,2,-nloctyp:nloctyp))
1559 allocate(e0new(3,-nloctyp:nloctyp))
1560 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1561 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1562 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1563 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1564 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1565 allocate(e0newtor(3,-nloctyp:nloctyp))
1566 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1568 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1569 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1570 itype2loc(ntyp1)=nloctyp
1571 iloctyp(nloctyp)=ntyp1
1573 itype2loc(-i)=-itype2loc(i)
1576 allocate(iloctyp(-nloctyp:nloctyp))
1577 allocate(itype2loc(-ntyp1:ntyp1))
1584 iloctyp(-i)=-iloctyp(i)
1586 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1587 !c write (iout,*) "nloctyp",nloctyp,
1588 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1589 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1590 !c write (iout,*) "nloctyp",nloctyp,
1591 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1594 !c write (iout,*) "NEWCORR",i
1595 read (ifourier,*,end=115,err=115)
1598 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1601 !c write (iout,*) "NEWCORR BNEW1"
1602 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1605 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1608 !c write (iout,*) "NEWCORR BNEW2"
1609 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1611 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1612 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1614 !c write (iout,*) "NEWCORR CCNEW"
1615 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1617 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1618 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1620 !c write (iout,*) "NEWCORR DDNEW"
1621 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1625 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1629 !c write (iout,*) "NEWCORR EENEW1"
1630 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1632 read (ifourier,*,end=115,err=115) e0new(ii,i)
1634 !c write (iout,*) (e0new(ii,i),ii=1,3)
1636 !c write (iout,*) "NEWCORR EENEW"
1639 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1640 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1641 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1642 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1647 bnew1(ii,1,-i)= bnew1(ii,1,i)
1648 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1649 bnew2(ii,1,-i)= bnew2(ii,1,i)
1650 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1653 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1654 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1655 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1656 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1657 ccnew(ii,1,-i)=ccnew(ii,1,i)
1658 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1659 ddnew(ii,1,-i)=ddnew(ii,1,i)
1660 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1662 e0new(1,-i)= e0new(1,i)
1663 e0new(2,-i)=-e0new(2,i)
1664 e0new(3,-i)=-e0new(3,i)
1666 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1667 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1668 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1669 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1673 write (iout,'(a)') "Coefficients of the multibody terms"
1674 do i=-nloctyp+1,nloctyp-1
1675 write (iout,*) "Type: ",onelet(iloctyp(i))
1676 write (iout,*) "Coefficients of the expansion of B1"
1678 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1680 write (iout,*) "Coefficients of the expansion of B2"
1682 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1684 write (iout,*) "Coefficients of the expansion of C"
1685 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1686 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1687 write (iout,*) "Coefficients of the expansion of D"
1688 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1689 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1690 write (iout,*) "Coefficients of the expansion of E"
1691 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1694 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1699 IF (SPLIT_FOURIERTOR) THEN
1701 !c write (iout,*) "NEWCORR TOR",i
1702 read (ifourier,*,end=115,err=115)
1705 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1708 !c write (iout,*) "NEWCORR BNEW1 TOR"
1709 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1712 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1715 !c write (iout,*) "NEWCORR BNEW2 TOR"
1716 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1718 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1719 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1721 !c write (iout,*) "NEWCORR CCNEW TOR"
1722 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1724 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1725 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1727 !c write (iout,*) "NEWCORR DDNEW TOR"
1728 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1732 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1736 !c write (iout,*) "NEWCORR EENEW1 TOR"
1737 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1739 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1741 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1743 !c write (iout,*) "NEWCORR EENEW TOR"
1746 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1747 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1748 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1749 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1754 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1755 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1756 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1757 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1760 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1761 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1762 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1763 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1765 e0newtor(1,-i)= e0newtor(1,i)
1766 e0newtor(2,-i)=-e0newtor(2,i)
1767 e0newtor(3,-i)=-e0newtor(3,i)
1769 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1770 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1771 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1772 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1776 write (iout,'(a)') &
1777 "Single-body coefficients of the torsional potentials"
1778 do i=-nloctyp+1,nloctyp-1
1779 write (iout,*) "Type: ",onelet(iloctyp(i))
1780 write (iout,*) "Coefficients of the expansion of B1tor"
1782 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1783 j,(bnew1tor(k,j,i),k=1,3)
1785 write (iout,*) "Coefficients of the expansion of B2tor"
1787 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1788 j,(bnew2tor(k,j,i),k=1,3)
1790 write (iout,*) "Coefficients of the expansion of Ctor"
1791 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1792 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1793 write (iout,*) "Coefficients of the expansion of Dtor"
1794 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1795 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1796 write (iout,*) "Coefficients of the expansion of Etor"
1797 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1800 write (iout,'(1hE,2i1,2f10.5)') &
1801 j,k,(eenewtor(l,j,k,i),l=1,2)
1807 do i=-nloctyp+1,nloctyp-1
1810 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1815 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1819 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1820 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1821 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1822 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1825 ENDIF !SPLIT_FOURIER_TOR
1827 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1828 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1829 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1830 allocate(b(13,-nloctyp-1:nloctyp+1))
1832 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1834 read (ifourier,*,end=115,err=115)
1835 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1837 write (iout,*) 'Type ',onelet(iloctyp(i))
1838 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1852 !c B1tilde(1,i) = b(3)
1853 !c! B1tilde(2,i) =-b(5)
1854 !c! B1tilde(1,-i) =-b(3)
1855 !c! B1tilde(2,-i) =b(5)
1856 !c! b1tilde(1,i)=0.0d0
1857 !c b1tilde(2,i)=0.0d0
1862 !cc B1tilde(1,i) = b(3,i)
1863 !cc B1tilde(2,i) =-b(5,i)
1864 !c B1tilde(1,-i) =-b(3,i)
1865 !c B1tilde(2,-i) =b(5,i)
1866 !cc b1tilde(1,i)=0.0d0
1867 !cc b1tilde(2,i)=0.0d0
1868 !cc B2(1,i) = b(2,i)
1869 !cc B2(2,i) = b(4,i)
1871 !c B2(2,-i) =-b(4,i)
1875 CCold(1,1,i)= b(7,i)
1876 CCold(2,2,i)=-b(7,i)
1877 CCold(2,1,i)= b(9,i)
1878 CCold(1,2,i)= b(9,i)
1879 CCold(1,1,-i)= b(7,i)
1880 CCold(2,2,-i)=-b(7,i)
1881 CCold(2,1,-i)=-b(9,i)
1882 CCold(1,2,-i)=-b(9,i)
1887 !c Ctilde(1,1,i)= CCold(1,1,i)
1888 !c Ctilde(1,2,i)= CCold(1,2,i)
1889 !c Ctilde(2,1,i)=-CCold(2,1,i)
1890 !c Ctilde(2,2,i)=-CCold(2,2,i)
1895 !c Ctilde(1,1,i)= CCold(1,1,i)
1896 !c Ctilde(1,2,i)= CCold(1,2,i)
1897 !c Ctilde(2,1,i)=-CCold(2,1,i)
1898 !c Ctilde(2,2,i)=-CCold(2,2,i)
1900 !c Ctilde(1,1,i)=0.0d0
1901 !c Ctilde(1,2,i)=0.0d0
1902 !c Ctilde(2,1,i)=0.0d0
1903 !c Ctilde(2,2,i)=0.0d0
1904 DDold(1,1,i)= b(6,i)
1905 DDold(2,2,i)=-b(6,i)
1906 DDold(2,1,i)= b(8,i)
1907 DDold(1,2,i)= b(8,i)
1908 DDold(1,1,-i)= b(6,i)
1909 DDold(2,2,-i)=-b(6,i)
1910 DDold(2,1,-i)=-b(8,i)
1911 DDold(1,2,-i)=-b(8,i)
1916 !c Dtilde(1,1,i)= DD(1,1,i)
1917 !c Dtilde(1,2,i)= DD(1,2,i)
1918 !c Dtilde(2,1,i)=-DD(2,1,i)
1919 !c Dtilde(2,2,i)=-DD(2,2,i)
1921 !c Dtilde(1,1,i)=0.0d0
1922 !c Dtilde(1,2,i)=0.0d0
1923 !c Dtilde(2,1,i)=0.0d0
1924 !c Dtilde(2,2,i)=0.0d0
1925 EEold(1,1,i)= b(10,i)+b(11,i)
1926 EEold(2,2,i)=-b(10,i)+b(11,i)
1927 EEold(2,1,i)= b(12,i)-b(13,i)
1928 EEold(1,2,i)= b(12,i)+b(13,i)
1929 EEold(1,1,-i)= b(10,i)+b(11,i)
1930 EEold(2,2,-i)=-b(10,i)+b(11,i)
1931 EEold(2,1,-i)=-b(12,i)+b(13,i)
1932 EEold(1,2,-i)=-b(12,i)-b(13,i)
1933 write(iout,*) "TU DOCHODZE"
1939 !c ee(2,1,i)=ee(1,2,i)
1944 "Coefficients of the cumulants (independent of valence angles)"
1945 do i=-nloctyp+1,nloctyp-1
1946 write (iout,*) 'Type ',onelet(iloctyp(i))
1948 write(iout,'(2f10.5)') B(3,i),B(5,i)
1950 write(iout,'(2f10.5)') B(2,i),B(4,i)
1953 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1957 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1961 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1970 ! Read torsional parameters in old format
1972 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1974 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1975 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1976 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1978 !el from energy module--------
1979 allocate(v1(nterm_old,ntortyp,ntortyp))
1980 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1981 !el---------------------------
1986 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1992 write (iout,'(/a/)') 'Torsional constants:'
1995 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1996 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2002 ! Read torsional parameters
2004 IF (TOR_MODE.eq.0) THEN
2005 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2006 read (itorp,*,end=113,err=113) ntortyp
2007 !el from energy module---------
2008 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2009 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2011 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2012 allocate(vlor2(maxlor,ntortyp,ntortyp))
2013 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2014 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2016 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2017 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2018 !el---------------------------
2021 !el---------------------------
2023 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2025 itortyp(i)=-itortyp(-i)
2027 itortyp(ntyp1)=ntortyp
2028 itortyp(-ntyp1)=-ntortyp
2030 write (iout,*) 'ntortyp',ntortyp
2032 do j=-ntortyp+1,ntortyp-1
2033 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2035 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2036 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2039 do k=1,nterm(i,j,iblock)
2040 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2042 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2043 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2044 v0ij=v0ij+si*v1(k,i,j,iblock)
2046 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2047 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2048 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2050 do k=1,nlor(i,j,iblock)
2051 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2052 vlor2(k,i,j),vlor3(k,i,j)
2053 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2056 v0(-i,-j,iblock)=v0ij
2062 write (iout,'(/a/)') 'Torsional constants:'
2064 do i=-ntortyp,ntortyp
2065 do j=-ntortyp,ntortyp
2066 write (iout,*) 'ityp',i,' jtyp',j
2067 write (iout,*) 'Fourier constants'
2068 do k=1,nterm(i,j,iblock)
2069 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2072 write (iout,*) 'Lorenz constants'
2073 do k=1,nlor(i,j,iblock)
2074 write (iout,'(3(1pe15.5))') &
2075 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2081 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2083 ! 6/23/01 Read parameters for double torsionals
2085 !el from energy module------------
2086 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2087 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2088 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2089 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2090 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2091 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2092 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2093 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2094 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2095 !---------------------------------
2099 do j=-ntortyp+1,ntortyp-1
2100 do k=-ntortyp+1,ntortyp-1
2101 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2102 ! write (iout,*) "OK onelett",
2105 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2106 .or. t3.ne.toronelet(k)) then
2107 write (iout,*) "Error in double torsional parameter file",&
2110 call MPI_Finalize(Ierror)
2112 stop "Error in double torsional parameter file"
2114 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2115 ntermd_2(i,j,k,iblock)
2116 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2117 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2118 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2119 ntermd_1(i,j,k,iblock))
2120 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2121 ntermd_1(i,j,k,iblock))
2122 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2123 ntermd_1(i,j,k,iblock))
2124 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2125 ntermd_1(i,j,k,iblock))
2126 ! Martix of D parameters for one dimesional foureir series
2127 do l=1,ntermd_1(i,j,k,iblock)
2128 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2129 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2130 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2131 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2132 ! write(iout,*) "whcodze" ,
2133 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2135 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2136 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2137 v2s(m,l,i,j,k,iblock),&
2138 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2139 ! Martix of D parameters for two dimesional fourier series
2140 do l=1,ntermd_2(i,j,k,iblock)
2142 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2143 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2144 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2145 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2154 write (iout,*) 'Constants for double torsionals'
2157 do j=-ntortyp+1,ntortyp-1
2158 do k=-ntortyp+1,ntortyp-1
2159 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2160 ' nsingle',ntermd_1(i,j,k,iblock),&
2161 ' ndouble',ntermd_2(i,j,k,iblock)
2163 write (iout,*) 'Single angles:'
2164 do l=1,ntermd_1(i,j,k,iblock)
2165 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2166 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2167 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2168 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2171 write (iout,*) 'Pairs of angles:'
2172 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2173 do l=1,ntermd_2(i,j,k,iblock)
2174 write (iout,'(i5,20f10.5)') &
2175 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2178 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2179 do l=1,ntermd_2(i,j,k,iblock)
2180 write (iout,'(i5,20f10.5)') &
2181 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2182 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2192 itype2loc(i)=itortyp(i)
2196 ELSE IF (TOR_MODE.eq.1) THEN
2198 !C read valence-torsional parameters
2199 read (itorp,*,end=121,err=121) ntortyp
2201 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2203 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2204 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2207 itype2loc(i)=itortyp(i)
2211 itortyp(i)=-itortyp(-i)
2213 do i=-ntortyp+1,ntortyp-1
2214 do j=-ntortyp+1,ntortyp-1
2215 !C first we read the cos and sin gamma parameters
2216 read (itorp,'(13x,a)',end=121,err=121) string
2217 write (iout,*) i,j,string
2218 read (itorp,*,end=121,err=121) &
2219 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2220 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2221 do k=1,nterm_kcc(j,i)
2222 do l=1,nterm_kcc_Tb(j,i)
2223 do ll=1,nterm_kcc_Tb(j,i)
2224 read (itorp,*,end=121,err=121) ii,jj,kk, &
2225 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2233 !c AL 4/8/16: Calculate coefficients from one-body parameters
2235 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2236 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2237 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2238 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2239 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2242 print *,i,itortyp(i)
2243 itortyp(i)=itype2loc(i)
2246 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2248 do i=-ntortyp+1,ntortyp-1
2249 do j=-ntortyp+1,ntortyp-1
2252 do k=1,nterm_kcc_Tb(j,i)
2253 do l=1,nterm_kcc_Tb(j,i)
2254 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2255 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2256 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2257 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2260 do k=1,nterm_kcc_Tb(j,i)
2261 do l=1,nterm_kcc_Tb(j,i)
2263 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2264 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2265 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2266 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2268 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2269 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2270 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2271 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2275 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2279 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2284 if (tor_mode.gt.0 .and. lprint) then
2285 !c Print valence-torsional parameters
2286 write (iout,'(a)') &
2287 "Parameters of the valence-torsional potentials"
2288 do i=-ntortyp+1,ntortyp-1
2289 do j=-ntortyp+1,ntortyp-1
2290 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2291 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2292 do k=1,nterm_kcc(j,i)
2293 do l=1,nterm_kcc_Tb(j,i)
2294 do ll=1,nterm_kcc_Tb(j,i)
2295 write (iout,'(3i5,2f15.4)')&
2296 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2304 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2305 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2306 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2307 !el from energy module---------
2308 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2309 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2311 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2312 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2313 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2314 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2316 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2317 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2318 !el---------------------------
2321 !el--------------------
2322 read (itorp_nucl,*,end=113,err=113) &
2323 (itortyp_nucl(i),i=1,ntyp_molec(2))
2324 ! print *,itortyp_nucl(:)
2325 !c write (iout,*) 'ntortyp',ntortyp
2328 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2329 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2332 do k=1,nterm_nucl(i,j)
2333 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2334 v0ij=v0ij+si*v1_nucl(k,i,j)
2337 do k=1,nlor_nucl(i,j)
2338 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2339 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2340 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2346 ! Read of Side-chain backbone correlation parameters
2347 ! Modified 11 May 2012 by Adasko
2350 read (isccor,*,end=119,err=119) nsccortyp
2352 !el from module energy-------------
2353 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2354 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2355 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2356 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2357 !-----------------------------------
2359 !el from module energy-------------
2360 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2362 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2364 isccortyp(i)=-isccortyp(-i)
2366 iscprol=isccortyp(20)
2367 ! write (iout,*) 'ntortyp',ntortyp
2369 !c maxinter is maximum interaction sites
2370 !el from module energy---------
2371 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2372 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2373 -nsccortyp:nsccortyp))
2374 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2375 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2376 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2377 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2378 !-----------------------------------
2382 read (isccor,*,end=119,err=119) &
2383 nterm_sccor(i,j),nlor_sccor(i,j)
2389 nterm_sccor(-i,j)=nterm_sccor(i,j)
2390 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2391 nterm_sccor(i,-j)=nterm_sccor(i,j)
2392 do k=1,nterm_sccor(i,j)
2393 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2395 if (j.eq.iscprol) then
2396 if (i.eq.isccortyp(10)) then
2397 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2398 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2400 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2401 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2402 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2403 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2404 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2405 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2406 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2407 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2410 if (i.eq.isccortyp(10)) then
2411 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2412 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2414 if (j.eq.isccortyp(10)) then
2415 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2416 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2418 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2419 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2420 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2421 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2422 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2423 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2427 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2428 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2429 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2430 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2433 do k=1,nlor_sccor(i,j)
2434 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2435 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2436 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2437 (1+vlor3sccor(k,i,j)**2)
2439 v0sccor(l,i,j)=v0ijsccor
2440 v0sccor(l,-i,j)=v0ijsccor1
2441 v0sccor(l,i,-j)=v0ijsccor2
2442 v0sccor(l,-i,-j)=v0ijsccor3
2448 !el from module energy-------------
2449 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2451 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2452 ! write (iout,*) 'ntortyp',ntortyp
2454 !c maxinter is maximum interaction sites
2455 !el from module energy---------
2456 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2457 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2458 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2459 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2460 !-----------------------------------
2464 read (isccor,*,end=119,err=119) &
2465 nterm_sccor(i,j),nlor_sccor(i,j)
2469 do k=1,nterm_sccor(i,j)
2470 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2472 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2475 do k=1,nlor_sccor(i,j)
2476 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2477 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2478 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2479 (1+vlor3sccor(k,i,j)**2)
2481 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2490 write (iout,'(/a/)') 'Torsional constants:'
2493 write (iout,*) 'ityp',i,' jtyp',j
2494 write (iout,*) 'Fourier constants'
2495 do k=1,nterm_sccor(i,j)
2496 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2498 write (iout,*) 'Lorenz constants'
2499 do k=1,nlor_sccor(i,j)
2500 write (iout,'(3(1pe15.5))') &
2501 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2508 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2509 ! interaction energy of the Gly, Ala, and Pro prototypes.
2512 ! Read electrostatic-interaction parameters
2517 write (iout,'(/a)') 'Electrostatic interaction constants:'
2518 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2519 'IT','JT','APP','BPP','AEL6','AEL3'
2521 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2522 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2523 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2524 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2529 app (i,j)=epp(i,j)*rri*rri
2530 bpp (i,j)=-2.0D0*epp(i,j)*rri
2531 ael6(i,j)=elpp6(i,j)*4.2D0**6
2532 ael3(i,j)=elpp3(i,j)*4.2D0**3
2534 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2540 ! Read side-chain interaction parameters.
2542 !el from module energy - COMMON.INTERACT-------
2543 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2544 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2545 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2547 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2548 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2549 allocate(epslip(ntyp,ntyp))
2557 !--------------------------------
2559 read (isidep,*,end=117,err=117) ipot,expon
2560 if (ipot.lt.1 .or. ipot.gt.5) then
2561 write (iout,'(2a)') 'Error while reading SC interaction',&
2562 'potential file - unknown potential type.'
2564 call MPI_Finalize(Ierror)
2570 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2571 ', exponents are ',expon,2*expon
2572 ! goto (10,20,30,30,40) ipot
2574 !----------------------- LJ potential ---------------------------------
2576 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2577 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2578 (sigma0(i),i=1,ntyp)
2580 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2581 write (iout,'(a/)') 'The epsilon array:'
2582 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2583 write (iout,'(/a)') 'One-body parameters:'
2584 write (iout,'(a,4x,a)') 'residue','sigma'
2585 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2588 !----------------------- LJK potential --------------------------------
2590 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2591 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2592 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2594 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2595 write (iout,'(a/)') 'The epsilon array:'
2596 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2597 write (iout,'(/a)') 'One-body parameters:'
2598 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2599 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2603 !---------------------- GB or BP potential -----------------------------
2606 ! print *,"I AM in SCELE",scelemode
2607 if (scelemode.eq.0) then
2609 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2611 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2612 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2613 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2614 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2616 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2619 ! For the GB potential convert sigma'**2 into chi'
2622 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2626 write (iout,'(/a/)') 'Parameters of the BP potential:'
2627 write (iout,'(a/)') 'The epsilon array:'
2628 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2629 write (iout,'(/a)') 'One-body parameters:'
2630 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2632 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2633 chip(i),alp(i),i=1,ntyp)
2636 ! print *,ntyp,"NTYP"
2637 allocate(icharge(ntyp1))
2638 ! print *,ntyp,icharge(i)
2640 read (isidep,*) (icharge(i),i=1,ntyp)
2641 print *,ntyp,icharge(i)
2642 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2643 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2644 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2645 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2646 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2647 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2648 allocate(epsintab(ntyp,ntyp))
2649 allocate(dtail(2,ntyp,ntyp))
2650 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2651 allocate(wqdip(2,ntyp,ntyp))
2652 allocate(wstate(4,ntyp,ntyp))
2653 allocate(dhead(2,2,ntyp,ntyp))
2654 allocate(nstate(ntyp,ntyp))
2655 allocate(debaykap(ntyp,ntyp))
2657 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2658 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2662 ! write (*,*) "Im in ALAB", i, " ", j
2664 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2665 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2666 chis(i,j),chis(j,i), & !2 w tej linii
2667 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2668 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2669 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2670 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2671 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2672 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2673 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2674 IF ((LaTeX).and.(i.gt.24)) then
2675 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2676 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2677 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2678 chis(i,j),chis(j,i) !2 w tej linii
2680 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2685 IF ((LaTeX).and.(i.gt.24)) then
2686 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2687 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2688 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2689 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2690 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2691 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2692 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2699 sigma(i,j) = sigma(j,i)
2700 sigmap1(i,j)=sigmap1(j,i)
2701 sigmap2(i,j)=sigmap2(j,i)
2702 sigiso1(i,j)=sigiso1(j,i)
2703 sigiso2(i,j)=sigiso2(j,i)
2704 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2705 nstate(i,j) = nstate(j,i)
2706 dtail(1,i,j) = dtail(1,j,i)
2707 dtail(2,i,j) = dtail(2,j,i)
2709 alphasur(k,i,j) = alphasur(k,j,i)
2710 wstate(k,i,j) = wstate(k,j,i)
2711 alphiso(k,i,j) = alphiso(k,j,i)
2714 dhead(2,1,i,j) = dhead(1,1,j,i)
2715 dhead(2,2,i,j) = dhead(1,2,j,i)
2716 dhead(1,1,i,j) = dhead(2,1,j,i)
2717 dhead(1,2,i,j) = dhead(2,2,j,i)
2719 epshead(i,j) = epshead(j,i)
2720 sig0head(i,j) = sig0head(j,i)
2723 wqdip(k,i,j) = wqdip(k,j,i)
2726 wquad(i,j) = wquad(j,i)
2727 epsintab(i,j) = epsintab(j,i)
2728 debaykap(i,j)=debaykap(j,i)
2729 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2734 !--------------------- GBV potential -----------------------------------
2736 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2737 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2738 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2739 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2741 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2742 write (iout,'(a/)') 'The epsilon array:'
2743 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2744 write (iout,'(/a)') 'One-body parameters:'
2745 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2746 's||/s_|_^2',' chip ',' alph '
2747 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2748 sigii(i),chip(i),alp(i),i=1,ntyp)
2751 write(iout,*)"Wrong ipot"
2757 !-----------------------------------------------------------------------
2758 ! Calculate the "working" parameters of SC interactions.
2760 !el from module energy - COMMON.INTERACT-------
2761 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2762 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2763 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2764 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2765 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2766 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2768 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2774 if (scelemode.eq.0) then
2783 sc_aa_tube_par(:)=0.0d0
2784 sc_bb_tube_par(:)=0.0d0
2786 !--------------------------------
2791 epslip(i,j)=epslip(j,i)
2794 if (scelemode.eq.0) then
2797 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2798 sigma(j,i)=sigma(i,j)
2799 rs0(i,j)=dwa16*sigma(i,j)
2804 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2805 'Working parameters of the SC interactions:',&
2806 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2811 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2813 ! print *,"SIGMA ZLA?",sigma(i,j)
2821 sigeps=dsign(1.0D0,epsij)
2823 aa_aq(i,j)=epsij*rrij*rrij
2824 ! print *,"ADASKO",epsij,rrij,expon
2825 bb_aq(i,j)=-sigeps*epsij*rrij
2826 aa_aq(j,i)=aa_aq(i,j)
2827 bb_aq(j,i)=bb_aq(i,j)
2828 epsijlip=epslip(i,j)
2829 sigeps=dsign(1.0D0,epsijlip)
2830 epsijlip=dabs(epsijlip)
2831 aa_lip(i,j)=epsijlip*rrij*rrij
2832 bb_lip(i,j)=-sigeps*epsijlip*rrij
2833 aa_lip(j,i)=aa_lip(i,j)
2834 bb_lip(j,i)=bb_lip(i,j)
2835 !C write(iout,*) aa_lip
2836 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2837 sigt1sq=sigma0(i)**2
2838 sigt2sq=sigma0(j)**2
2841 ratsig1=sigt2sq/sigt1sq
2842 ratsig2=1.0D0/ratsig1
2843 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2844 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2845 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2849 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2850 sigmaii(i,j)=rsum_max
2851 sigmaii(j,i)=rsum_max
2853 ! sigmaii(i,j)=r0(i,j)
2854 ! sigmaii(j,i)=r0(i,j)
2856 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2857 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2858 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2859 augm(i,j)=epsij*r_augm**(2*expon)
2860 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2867 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2868 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2869 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2874 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2875 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2876 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2877 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2878 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2879 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2880 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2881 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2882 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2883 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2884 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2885 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2886 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2887 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2888 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2889 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2898 read (isidep_nucl,*) ipot_nucl
2899 ! print *,"TU?!",ipot_nucl
2900 if (ipot_nucl.eq.1) then
2901 do i=1,ntyp_molec(2)
2902 do j=i,ntyp_molec(2)
2903 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2904 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2908 do i=1,ntyp_molec(2)
2909 do j=i,ntyp_molec(2)
2910 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2911 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2912 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2916 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2917 do i=1,ntyp_molec(2)
2918 do j=i,ntyp_molec(2)
2919 rrij=sigma_nucl(i,j)
2923 epsij=4*eps_nucl(i,j)
2924 sigeps=dsign(1.0D0,epsij)
2926 aa_nucl(i,j)=epsij*rrij*rrij
2927 bb_nucl(i,j)=-sigeps*epsij*rrij
2928 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2929 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2930 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2931 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2932 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2933 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2936 aa_nucl(i,j)=aa_nucl(j,i)
2937 bb_nucl(i,j)=bb_nucl(j,i)
2938 ael3_nucl(i,j)=ael3_nucl(j,i)
2939 ael6_nucl(i,j)=ael6_nucl(j,i)
2940 ael63_nucl(i,j)=ael63_nucl(j,i)
2941 ael32_nucl(i,j)=ael32_nucl(j,i)
2942 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2943 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2944 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2945 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2946 eps_nucl(i,j)=eps_nucl(j,i)
2947 sigma_nucl(i,j)=sigma_nucl(j,i)
2948 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2952 write(iout,*) "tube param"
2953 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2954 ccavtubpep,dcavtubpep,tubetranenepep
2955 sigmapeptube=sigmapeptube**6
2956 sigeps=dsign(1.0D0,epspeptube)
2957 epspeptube=dabs(epspeptube)
2958 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2959 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2960 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2962 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2963 ccavtub(i),dcavtub(i),tubetranene(i)
2964 sigmasctube=sigmasctube**6
2965 sigeps=dsign(1.0D0,epssctube)
2966 epssctube=dabs(epssctube)
2967 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2968 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2969 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2971 !-----------------READING SC BASE POTENTIALS-----------------------------
2972 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2973 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2974 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2975 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2976 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2977 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2978 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2979 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2980 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2981 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2982 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2983 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2984 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2985 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2986 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2987 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2990 do i=1,ntyp_molec(1)
2991 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2992 write (*,*) "Im in ", i, " ", j
2993 read(isidep_scbase,*) &
2994 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2995 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2996 write(*,*) "eps",eps_scbase(i,j)
2997 read(isidep_scbase,*) &
2998 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2999 chis_scbase(i,j,1),chis_scbase(i,j,2)
3000 read(isidep_scbase,*) &
3001 dhead_scbasei(i,j), &
3002 dhead_scbasej(i,j), &
3003 rborn_scbasei(i,j),rborn_scbasej(i,j)
3004 read(isidep_scbase,*) &
3005 (wdipdip_scbase(k,i,j),k=1,3), &
3006 (wqdip_scbase(k,i,j),k=1,2)
3007 read(isidep_scbase,*) &
3008 alphapol_scbase(i,j), &
3009 epsintab_scbase(i,j)
3012 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3013 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3015 do i=1,ntyp_molec(1)
3016 do j=1,ntyp_molec(2)-1
3017 epsij=eps_scbase(i,j)
3018 rrij=sigma_scbase(i,j)
3023 sigeps=dsign(1.0D0,epsij)
3025 aa_scbase(i,j)=epsij*rrij*rrij
3026 bb_scbase(i,j)=-sigeps*epsij*rrij
3029 !-----------------READING PEP BASE POTENTIALS-------------------
3030 allocate(eps_pepbase(ntyp_molec(2)))
3031 allocate(sigma_pepbase(ntyp_molec(2)))
3032 allocate(chi_pepbase(ntyp_molec(2),2))
3033 allocate(chipp_pepbase(ntyp_molec(2),2))
3034 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3035 allocate(sigmap1_pepbase(ntyp_molec(2)))
3036 allocate(sigmap2_pepbase(ntyp_molec(2)))
3037 allocate(chis_pepbase(ntyp_molec(2),2))
3038 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3041 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3042 write (*,*) "Im in ", i, " ", j
3043 read(isidep_pepbase,*) &
3044 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3045 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3046 write(*,*) "eps",eps_pepbase(j)
3047 read(isidep_pepbase,*) &
3048 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3049 chis_pepbase(j,1),chis_pepbase(j,2)
3050 read(isidep_pepbase,*) &
3051 (wdipdip_pepbase(k,j),k=1,3)
3053 allocate(aa_pepbase(ntyp_molec(2)))
3054 allocate(bb_pepbase(ntyp_molec(2)))
3056 do j=1,ntyp_molec(2)-1
3057 epsij=eps_pepbase(j)
3058 rrij=sigma_pepbase(j)
3063 sigeps=dsign(1.0D0,epsij)
3065 aa_pepbase(j)=epsij*rrij*rrij
3066 bb_pepbase(j)=-sigeps*epsij*rrij
3068 !--------------READING SC PHOSPHATE-------------------------------------
3069 allocate(eps_scpho(ntyp_molec(1)))
3070 allocate(sigma_scpho(ntyp_molec(1)))
3071 allocate(chi_scpho(ntyp_molec(1),2))
3072 allocate(chipp_scpho(ntyp_molec(1),2))
3073 allocate(alphasur_scpho(4,ntyp_molec(1)))
3074 allocate(sigmap1_scpho(ntyp_molec(1)))
3075 allocate(sigmap2_scpho(ntyp_molec(1)))
3076 allocate(chis_scpho(ntyp_molec(1),2))
3077 allocate(wqq_scpho(ntyp_molec(1)))
3078 allocate(wqdip_scpho(2,ntyp_molec(1)))
3079 allocate(alphapol_scpho(ntyp_molec(1)))
3080 allocate(epsintab_scpho(ntyp_molec(1)))
3081 allocate(dhead_scphoi(ntyp_molec(1)))
3082 allocate(rborn_scphoi(ntyp_molec(1)))
3083 allocate(rborn_scphoj(ntyp_molec(1)))
3084 allocate(alphi_scpho(ntyp_molec(1)))
3088 do j=1,ntyp_molec(1) ! without U then we will take T for U
3089 write (*,*) "Im in scpho ", i, " ", j
3090 read(isidep_scpho,*) &
3091 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3092 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3093 write(*,*) "eps",eps_scpho(j)
3094 read(isidep_scpho,*) &
3095 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3096 chis_scpho(j,1),chis_scpho(j,2)
3097 read(isidep_scpho,*) &
3098 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3099 read(isidep_scpho,*) &
3100 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3104 allocate(aa_scpho(ntyp_molec(1)))
3105 allocate(bb_scpho(ntyp_molec(1)))
3107 do j=1,ntyp_molec(1)
3114 sigeps=dsign(1.0D0,epsij)
3116 aa_scpho(j)=epsij*rrij*rrij
3117 bb_scpho(j)=-sigeps*epsij*rrij
3121 read(isidep_peppho,*) &
3122 eps_peppho,sigma_peppho
3123 read(isidep_peppho,*) &
3124 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3125 read(isidep_peppho,*) &
3126 (wqdip_peppho(k),k=1,2)
3134 sigeps=dsign(1.0D0,epsij)
3136 aa_peppho=epsij*rrij*rrij
3137 bb_peppho=-sigeps*epsij*rrij
3140 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3145 ! Define the SC-p interaction constants (hard-coded; old style)
3148 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3150 ! aad(i,1)=0.3D0*4.0D0**12
3151 ! Following line for constants currently implemented
3152 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3153 aad(i,1)=1.5D0*4.0D0**12
3154 ! aad(i,1)=0.17D0*5.6D0**12
3156 ! "Soft" SC-p repulsion
3158 ! Following line for constants currently implemented
3159 ! aad(i,1)=0.3D0*4.0D0**6
3160 ! "Hard" SC-p repulsion
3161 bad(i,1)=3.0D0*4.0D0**6
3162 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3171 ! 8/9/01 Read the SC-p interaction constants from file
3174 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3177 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3178 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3179 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3180 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3184 write (iout,*) "Parameters of SC-p interactions:"
3186 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3187 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3192 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3194 do i=1,ntyp_molec(2)
3195 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3197 do i=1,ntyp_molec(2)
3198 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3199 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3201 r0pp=1.12246204830937298142*5.16158
3207 ! Define the constants of the disulfide bridge
3211 ! Old arbitrary potential - commented out.
3216 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3217 ! energy surface of diethyl disulfide.
3218 ! A. Liwo and U. Kozlowska, 11/24/03
3235 write (iout,'(/a)') "Disulfide bridge parameters:"
3236 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3237 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3238 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3239 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3242 if (shield_mode.gt.0) then
3243 pi=4.0D0*datan(1.0D0)
3244 !C VSolvSphere the volume of solving sphere
3246 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3247 !C there will be no distinction between proline peptide group and normal peptide
3248 !C group in case of shielding parameters
3249 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3250 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3251 write (iout,*) VSolvSphere,VSolvSphere_div
3252 !C long axis of side chain
3254 long_r_sidechain(i)=vbldsc0(1,i)
3255 ! if (scelemode.eq.0) then
3256 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3257 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3259 ! short_r_sidechain(i)=sigma(i,i)
3261 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3268 111 write (iout,*) "Error reading bending energy parameters."
3270 112 write (iout,*) "Error reading rotamer energy parameters."
3272 113 write (iout,*) "Error reading torsional energy parameters."
3274 114 write (iout,*) "Error reading double torsional energy parameters."
3276 115 write (iout,*) &
3277 "Error reading cumulant (multibody energy) parameters."
3279 116 write (iout,*) "Error reading electrostatic energy parameters."
3281 117 write (iout,*) "Error reading side chain interaction parameters."
3283 118 write (iout,*) "Error reading SCp interaction parameters."
3285 119 write (iout,*) "Error reading SCCOR parameters"
3287 121 write (iout,*) "Error in Czybyshev parameters"
3290 call MPI_Finalize(Ierror)
3294 end subroutine parmread
3296 !-----------------------------------------------------------------------------
3298 !-----------------------------------------------------------------------------
3299 subroutine printmat(ldim,m,n,iout,key,a)
3302 character(len=3),dimension(n) :: key
3303 real(kind=8),dimension(ldim,n) :: a
3305 integer :: i,j,k,m,iout,nlim
3309 write (iout,1000) (key(k),k=i,nlim)
3311 1000 format (/5x,8(6x,a3))
3312 1020 format (/80(1h-)/)
3314 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3317 1010 format (a3,2x,8(f9.4))
3319 end subroutine printmat
3320 !-----------------------------------------------------------------------------
3322 !-----------------------------------------------------------------------------
3324 ! Read the PDB file and convert the peptide geometry into virtual-chain
3327 use energy_data, only: itype,istype
3331 ! use control, only: rescode,sugarcode
3332 ! implicit real*8 (a-h,o-z)
3333 ! include 'DIMENSIONS'
3334 ! include 'COMMON.LOCAL'
3335 ! include 'COMMON.VAR'
3336 ! include 'COMMON.CHAIN'
3337 ! include 'COMMON.INTERACT'
3338 ! include 'COMMON.IOUNITS'
3339 ! include 'COMMON.GEO'
3340 ! include 'COMMON.NAMES'
3341 ! include 'COMMON.CONTROL'
3342 ! include 'COMMON.DISTFIT'
3343 ! include 'COMMON.SETUP'
3344 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3346 logical :: lprn=.true.,fail
3347 real(kind=8),dimension(3) :: e1,e2,e3
3348 real(kind=8) :: dcj,efree_temp
3349 character(len=3) :: seq,res,res2
3350 character(len=5) :: atom
3351 character(len=80) :: card
3352 real(kind=8),dimension(3,20) :: sccor
3353 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3354 integer :: isugar,molecprev,firstion
3355 character*1 :: sugar
3357 real(kind=8),dimension(3) :: ccc
3359 integer,dimension(2,maxres/3) :: hfrag_alloc
3360 integer,dimension(4,maxres/3) :: bfrag_alloc
3361 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3362 real(kind=8),dimension(:,:), allocatable :: c_temporary
3363 integer,dimension(:,:) , allocatable :: itype_temporary
3364 integer,dimension(:),allocatable :: istype_temp
3371 ! write (2,*) "UNRES_PDB",unres_pdb
3391 !-----------------------------
3392 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3393 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3394 if(.not. allocated(istype)) allocate(istype(maxres))
3396 read (ipdbin,'(a80)',end=10) card
3397 write (iout,'(a)') card
3398 if (card(:5).eq.'HELIX') then
3401 read(card(22:25),*) hfrag(1,nhfrag)
3402 read(card(34:37),*) hfrag(2,nhfrag)
3404 if (card(:5).eq.'SHEET') then
3407 read(card(24:26),*) bfrag(1,nbfrag)
3408 read(card(35:37),*) bfrag(2,nbfrag)
3409 !rc----------------------------------------
3410 !rc to be corrected !!!
3411 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3412 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3413 !rc----------------------------------------
3415 if (card(:3).eq.'END') then
3417 else if (card(:3).eq.'TER') then
3422 itype(ires_old,molecule)=ntyp1_molec(molecule)
3423 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3424 nres_molec(molecule)=nres_molec(molecule)+2
3426 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3429 dc(j,ires)=sccor(j,iii)
3432 call sccenter(ires,iii,sccor)
3438 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3439 ! Fish out the ATOM cards.
3440 ! write(iout,*) 'card',card(1:20)
3441 ! print *,"ATU ",card(1:6), CARD(3:6)
3443 if (index(card(1:4),'ATOM').gt.0) then
3444 read (card(12:16),*) atom
3445 ! write (iout,*) "! ",atom," !",ires
3446 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3447 read (card(23:26),*) ires
3448 read (card(18:20),'(a3)') res
3449 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3450 ! & " ires_old",ires_old
3451 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3452 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3453 if (ires-ishift+ishift1.ne.ires_old) then
3454 ! Calculate the CM of the preceding residue.
3455 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3457 ! write (iout,*) "Calculating sidechain center iii",iii
3460 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3463 call sccenter(ires_old,iii,sccor)
3467 ! Start new residue.
3468 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3471 else if (ibeg.eq.1) then
3472 write (iout,*) "BEG ires",ires
3474 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3477 nres_molec(molecule)=nres_molec(molecule)+1
3479 ires=ires-ishift+ishift1
3481 ! write (iout,*) "ishift",ishift," ires",ires,&
3482 ! " ires_old",ires_old
3484 else if (ibeg.eq.2) then
3486 ishift=-ires_old+ires-1 !!!!!
3487 ishift1=ishift1-1 !!!!!
3488 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3489 ires=ires-ishift+ishift1
3490 ! print *,ires,ishift,ishift1
3494 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3495 ires=ires-ishift+ishift1
3498 ! print *,'atom',ires,atom
3499 if (res.eq.'ACE' .or. res.eq.'NHE') then
3502 if (atom.eq.'CA '.or.atom.eq.'N ') then
3504 itype(ires,molecule)=rescode(ires,res,0,molecule)
3506 ! nres_molec(molecule)=nres_molec(molecule)+1
3510 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3511 ! nres_molec(molecule)=nres_molec(molecule)+1
3512 read (card(19:19),'(a1)') sugar
3513 isugar=sugarcode(sugar,ires)
3514 ! if (ibeg.eq.1) then
3518 ! print *,"ires=",ires,istype(ires)
3524 ires=ires-ishift+ishift1
3526 ! write (iout,*) "ires_old",ires_old," ires",ires
3527 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3530 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3531 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3532 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3533 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3534 ! print *,ires,ishift,ishift1
3535 ! write (iout,*) "backbone ",atom
3537 write (iout,'(2i3,2x,a,3f8.3)') &
3538 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3541 nres_molec(molecule)=nres_molec(molecule)+1
3543 sccor(j,iii)=c(j,ires)
3545 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3546 atom.eq."C2'" .or. atom.eq."C3'" &
3547 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3548 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3549 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3550 ! print *,ires,ishift,ishift1
3554 ! sccor(j,iii)=c(j,ires)
3557 c(j,ires)=c(j,ires)+ccc(j)/5.0
3559 print *,counter,molecule
3560 if (counter.eq.5) then
3562 nres_molec(molecule)=nres_molec(molecule)+1
3565 ! sccor(j,iii)=c(j,ires)
3569 ! print *, "ATOM",atom(1:3)
3570 else if (atom.eq."C5'") then
3571 read (card(19:19),'(a1)') sugar
3572 isugar=sugarcode(sugar,ires)
3577 ! print *,ires,istype(ires)
3581 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3582 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3583 nres_molec(molecule)=nres_molec(molecule)+1
3584 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3588 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3590 else if ((atom.eq."C1'").and.unres_pdb) then
3592 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3593 ! write (*,*) card(23:27),ires,itype(ires,1)
3594 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3595 atom.ne.'N' .and. atom.ne.'C' .and. &
3596 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3597 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3598 .and. atom.ne.'P '.and. &
3599 atom(1:1).ne.'H' .and. &
3600 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3602 ! write (iout,*) "sidechain ",atom
3603 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3604 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3605 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3607 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3610 ! print *,"IONS",ions,card(1:6)
3611 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3612 if (firstion.eq.0) then
3616 dc(j,ires)=sccor(j,iii)
3619 call sccenter(ires,iii,sccor)
3622 read (card(12:16),*) atom
3623 ! print *,"HETATOM", atom
3624 read (card(18:20),'(a3)') res
3625 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3626 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3627 .or.(atom(1:2).eq.'K ')) &
3630 if (molecule.ne.5) molecprev=molecule
3632 nres_molec(molecule)=nres_molec(molecule)+1
3633 print *,"HERE",nres_molec(molecule)
3635 itype(ires,molecule)=rescode(ires,res,0,molecule)
3636 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3640 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3641 if (ires.eq.0) return
3642 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3645 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3647 nres_molec(molecule)=nres_molec(molecule)-2
3648 print *,'I have',nres, nres_molec(:)
3650 do k=1,4 ! ions are without dummy
3651 if (nres_molec(k).eq.0) cycle
3653 ! write (iout,*) i,itype(i,1)
3654 ! if (itype(i,1).eq.ntyp1) then
3655 ! write (iout,*) "dummy",i,itype(i,1)
3657 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3658 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3662 if (itype(i,k).eq.ntyp1_molec(k)) then
3663 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3664 if (itype(i+2,k).eq.0) then
3665 ! print *,"masz sieczke"
3667 if (itype(i+2,j).ne.0) then
3669 itype(i+1,j)=ntyp1_molec(j)
3670 nres_molec(k)=nres_molec(k)-1
3671 nres_molec(j)=nres_molec(j)+1
3677 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3678 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3679 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3680 ! if (unres_pdb) then
3681 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3682 ! print *,i,'tu dochodze'
3683 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3691 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3695 dcj=(c(j,i-2)-c(j,i-3))/2.0
3696 if (dcj.eq.0) dcj=1.23591524223
3701 else !itype(i+1,1).eq.ntyp1
3702 ! if (unres_pdb) then
3703 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3704 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3711 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3712 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3716 dcj=(c(j,i+3)-c(j,i+2))/2.0
3717 if (dcj.eq.0) dcj=1.23591524223
3722 endif !itype(i+1,1).eq.ntyp1
3723 endif !itype.eq.ntyp1
3727 ! Calculate the CM of the last side chain.
3731 dc(j,ires)=sccor(j,iii)
3734 call sccenter(ires,iii,sccor)
3740 ! print *,"molecule",molecule
3741 if ((itype(nres,1).ne.10)) then
3743 if (molecule.eq.5) molecule=molecprev
3744 itype(nres,molecule)=ntyp1_molec(molecule)
3745 nres_molec(molecule)=nres_molec(molecule)+1
3746 ! if (unres_pdb) then
3747 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3748 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3755 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3759 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3760 c(j,nres)=c(j,nres-1)+dcj
3761 c(j,2*nres)=c(j,nres)
3765 ! print *,'I have',nres, nres_molec(:)
3767 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3769 if (nres.ne.nres0) then
3770 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3772 stop "Error nres value in WHAM input"
3775 !---------------------------------
3776 !el reallocate tables
3779 ! hfrag_alloc(j,i)=hfrag(j,i)
3782 ! bfrag_alloc(j,i)=bfrag(j,i)
3788 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3789 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3790 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3791 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3795 ! hfrag(j,i)=hfrag_alloc(j,i)
3800 ! bfrag(j,i)=bfrag_alloc(j,i)
3803 !el end reallocate tables
3804 !---------------------------------
3812 c(j,2*nres)=c(j,nres)
3815 if (itype(1,1).eq.ntyp1) then
3819 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3820 call refsys(2,3,4,e1,e2,e3,fail)
3827 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3828 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
3832 dcj=(c(j,4)-c(j,3))/2.0
3838 ! First lets assign correct dummy to correct type of chain
3840 if (itype(1,1).eq.ntyp1) then
3841 if (itype(2,1).eq.0) then
3843 if (itype(2,j).ne.0) then
3845 itype(1,j)=ntyp1_molec(j)
3846 nres_molec(1)=nres_molec(1)-1
3847 nres_molec(j)=nres_molec(j)+1
3854 print *,'I have',nres, nres_molec(:)
3856 ! Copy the coordinates to reference coordinates
3862 ! Calculate internal coordinates.
3864 write (iout,'(/a)') &
3865 "Cartesian coordinates of the reference structure"
3866 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3867 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3869 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3870 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3871 (c(j,ires+nres),j=1,3)
3874 ! znamy już nres wiec mozna alokowac tablice
3875 ! Calculate internal coordinates.
3876 if(me.eq.king.or..not.out1file)then
3877 write (iout,'(a)') &
3878 "Backbone and SC coordinates as read from the PDB"
3880 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
3881 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3882 (c(j,nres+ires),j=1,3)
3885 ! NOW LETS ROCK! SORTING
3886 allocate(c_temporary(3,2*nres))
3887 allocate(itype_temporary(nres,5))
3888 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3889 if (.not.allocated(istype)) write(iout,*) &
3890 "SOMETHING WRONG WITH ISTYTPE"
3891 allocate(istype_temp(nres))
3892 itype_temporary(:,:)=0
3896 if (itype(i,k).ne.0) then
3898 c_temporary(j,seqalingbegin)=c(j,i)
3899 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3902 itype_temporary(seqalingbegin,k)=itype(i,k)
3903 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3904 istype_temp(seqalingbegin)=istype(i)
3905 molnum(seqalingbegin)=k
3906 seqalingbegin=seqalingbegin+1
3912 c(j,i)=c_temporary(j,i)
3917 itype(i,k)=itype_temporary(i,k)
3918 istype(i)=istype_temp(i)
3921 ! if (itype(1,1).eq.ntyp1) then
3924 ! if (unres_pdb) then
3925 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3926 ! call refsys(2,3,4,e1,e2,e3,fail)
3933 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3937 ! dcj=(c(j,4)-c(j,3))/2.0
3939 ! c(j,nres+1)=c(j,1)
3945 write (iout,'(/a)') &
3946 "Cartesian coordinates of the reference structure after sorting"
3947 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3948 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3950 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3951 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3952 (c(j,ires+nres),j=1,3)
3956 ! print *,seqalingbegin,nres
3957 if(.not.allocated(vbld)) then
3958 allocate(vbld(2*nres))
3963 if(.not.allocated(vbld_inv)) then
3964 allocate(vbld_inv(2*nres))
3970 if(.not.allocated(theta)) then
3971 allocate(theta(nres+2))
3975 if(.not.allocated(phi)) allocate(phi(nres+2))
3976 if(.not.allocated(alph)) allocate(alph(nres+2))
3977 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3978 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3979 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3980 if(.not.allocated(costtab)) allocate(costtab(nres))
3981 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3982 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3983 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3984 if(.not.allocated(xxref)) allocate(xxref(nres))
3985 if(.not.allocated(yyref)) allocate(yyref(nres))
3986 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3987 if(.not.allocated(dc_norm)) then
3988 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3989 allocate(dc_norm(3,0:2*nres+2))
3993 call int_from_cart(.true.,.false.)
3994 call sc_loc_geom(.false.)
3996 thetaref(i)=theta(i)
4006 dc(j,i)=c(j,i+1)-c(j,i)
4007 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4012 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4013 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4015 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4019 ! Copy the coordinates to reference coordinates
4020 ! Splits to single chain if occurs
4024 ! cref(j,i,cou)=c(j,i)
4028 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4029 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4030 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4031 !-----------------------------
4035 write (iout,*) "symetr", symetr
4038 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4040 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4043 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4048 cref(j,i,cou)=c(j,i)
4049 cref(j,i+nres,cou)=c(j,i+nres)
4051 chain_rep(j,lll,kkk)=c(j,i)
4052 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4056 write (iout,*) chain_length
4057 if (chain_length.eq.0) chain_length=nres
4059 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4060 chain_rep(j,chain_length+nres,symetr) &
4061 =chain_rep(j,chain_length+nres,1)
4064 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4066 ! do kkk=1,chain_length
4067 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4071 ! makes copy of chains
4072 write (iout,*) "symetr", symetr
4077 if (symetr.gt.1) then
4084 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4090 write (iout,*) i,icha
4091 do lll=1,chain_length
4093 if (cou.le.nres) then
4095 kupa=mod(lll,chain_length)
4096 iprzes=(kkk-1)*chain_length+lll
4097 if (kupa.eq.0) kupa=chain_length
4098 write (iout,*) "kupa", kupa
4099 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4100 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4107 !-koniec robienia kopii
4110 write (iout,*) "nowa struktura", nperm
4112 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4114 cref(3,i,kkk),cref(1,nres+i,kkk),&
4115 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4117 100 format (//' alpha-carbon coordinates ',&
4118 ' centroid coordinates'/ &
4119 ' ', 6X,'X',11X,'Y',11X,'Z', &
4120 10X,'X',11X,'Y',11X,'Z')
4121 110 format (a,'(',i5,')',6f12.5)
4127 bfrag(i,j)=bfrag(i,j)-ishift
4133 hfrag(i,j)=hfrag(i,j)-ishift
4139 end subroutine readpdb
4140 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4141 !-----------------------------------------------------------------------------
4143 !-----------------------------------------------------------------------------
4144 subroutine read_control
4158 use random, only: random_init
4159 ! implicit real*8 (a-h,o-z)
4160 ! include 'DIMENSIONS'
4162 use prng, only:prng_restart
4164 logical :: OKRandom!, prng_restart
4167 ! include 'COMMON.IOUNITS'
4168 ! include 'COMMON.TIME1'
4169 ! include 'COMMON.THREAD'
4170 ! include 'COMMON.SBRIDGE'
4171 ! include 'COMMON.CONTROL'
4172 ! include 'COMMON.MCM'
4173 ! include 'COMMON.MAP'
4174 ! include 'COMMON.HEADER'
4175 ! include 'COMMON.CSA'
4176 ! include 'COMMON.CHAIN'
4177 ! include 'COMMON.MUCA'
4178 ! include 'COMMON.MD'
4179 ! include 'COMMON.FFIELD'
4180 ! include 'COMMON.INTERACT'
4181 ! include 'COMMON.SETUP'
4182 !el integer :: KDIAG,ICORFL,IXDR
4183 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4184 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4185 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4186 ! character(len=80) :: ucase
4187 character(len=640) :: controlcard
4189 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4195 read (INP,'(a)') titel
4196 call card_concat(controlcard,.true.)
4197 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4198 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4199 call reada(controlcard,'SEED',seed,0.0D0)
4200 call random_init(seed)
4201 ! Set up the time limit (caution! The time must be input in minutes!)
4202 read_cart=index(controlcard,'READ_CART').gt.0
4203 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4204 call readi(controlcard,'SYM',symetr,1)
4205 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4206 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4207 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4208 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4209 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4210 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4211 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4212 call reada(controlcard,'DRMS',drms,0.1D0)
4213 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4214 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4215 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4216 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4217 write (iout,'(a,f10.1)')'DRMS = ',drms
4218 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4219 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4221 call readi(controlcard,'NZ_START',nz_start,0)
4222 call readi(controlcard,'NZ_END',nz_end,0)
4223 call readi(controlcard,'IZ_SC',iz_sc,0)
4224 timlim=60.0D0*timlim
4225 safety = 60.0d0*safety
4228 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4229 !C SHIELD keyword sets if the shielding effect of side-chains is used
4230 !C 0 denots no shielding is used all peptide are equally despite the
4231 !C solvent accesible area
4232 !C 1 the newly introduced function
4233 !C 2 reseved for further possible developement
4234 call readi(controlcard,'SHIELD',shield_mode,0)
4235 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4236 write(iout,*) "shield_mode",shield_mode
4237 call readi(controlcard,'TORMODE',tor_mode,0)
4238 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4239 write(iout,*) "torsional and valence angle mode",tor_mode
4241 !C Varibles set size of box
4242 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4243 protein=index(controlcard,"PROTEIN").gt.0
4244 ions=index(controlcard,"IONS").gt.0
4245 nucleic=index(controlcard,"NUCLEIC").gt.0
4246 write (iout,*) "with_theta_constr ",with_theta_constr
4247 AFMlog=(index(controlcard,'AFM'))
4248 selfguide=(index(controlcard,'SELFGUIDE'))
4249 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4250 call readi(controlcard,'GENCONSTR',genconstr,0)
4251 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4252 call reada(controlcard,'BOXY',boxysize,100.0d0)
4253 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4254 call readi(controlcard,'TUBEMOD',tubemode,0)
4255 print *,"SCELE",scelemode
4256 call readi(controlcard,"SCELEMODE",scelemode,0)
4257 print *,"SCELE",scelemode
4259 ! elemode = 0 is orignal UNRES electrostatics
4260 ! elemode = 1 is "Momo" potentials in progress
4261 ! elemode = 2 is in development EVALD
4264 write (iout,*) TUBEmode,"TUBEMODE"
4265 if (TUBEmode.gt.0) then
4266 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4267 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4268 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4269 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4270 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4271 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4272 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4273 buftubebot=bordtubebot+tubebufthick
4274 buftubetop=bordtubetop-tubebufthick
4277 ! CUTOFFF ON ELECTROSTATICS
4278 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
4279 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4280 write(iout,*) "R_CUT_ELE=",r_cut_ele
4281 ! Lipidic parameters
4282 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4283 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4284 if (lipthick.gt.0.0d0) then
4285 bordliptop=(boxzsize+lipthick)/2.0
4286 bordlipbot=bordliptop-lipthick
4287 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4288 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4289 buflipbot=bordlipbot+lipbufthick
4290 bufliptop=bordliptop-lipbufthick
4291 if ((lipbufthick*2.0d0).gt.lipthick) &
4292 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4293 endif !lipthick.gt.0
4294 write(iout,*) "bordliptop=",bordliptop
4295 write(iout,*) "bordlipbot=",bordlipbot
4296 write(iout,*) "bufliptop=",bufliptop
4297 write(iout,*) "buflipbot=",buflipbot
4298 write (iout,*) "SHIELD MODE",shield_mode
4300 !C-------------------------
4301 minim=(index(controlcard,'MINIMIZE').gt.0)
4302 dccart=(index(controlcard,'CART').gt.0)
4303 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4304 overlapsc=.not.overlapsc
4305 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4306 searchsc=.not.searchsc
4307 sideadd=(index(controlcard,'SIDEADD').gt.0)
4308 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4309 outpdb=(index(controlcard,'PDBOUT').gt.0)
4310 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4311 pdbref=(index(controlcard,'PDBREF').gt.0)
4312 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4313 indpdb=index(controlcard,'PDBSTART')
4314 extconf=(index(controlcard,'EXTCONF').gt.0)
4315 call readi(controlcard,'IPRINT',iprint,0)
4316 call readi(controlcard,'MAXGEN',maxgen,10000)
4317 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4318 call readi(controlcard,"KDIAG",kdiag,0)
4319 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4320 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4321 write (iout,*) "RESCALE_MODE",rescale_mode
4322 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4323 if (index(controlcard,'REGULAR').gt.0.0D0) then
4324 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4328 if (index(controlcard,'CHECKGRAD').gt.0) then
4330 if (index(controlcard,'CART').gt.0) then
4332 elseif (index(controlcard,'CARINT').gt.0) then
4337 elseif (index(controlcard,'THREAD').gt.0) then
4339 call readi(controlcard,'THREAD',nthread,0)
4340 if (nthread.gt.0) then
4341 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4344 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4345 stop 'Error termination in Read_Control.'
4347 else if (index(controlcard,'MCMA').gt.0) then
4349 else if (index(controlcard,'MCEE').gt.0) then
4351 else if (index(controlcard,'MULTCONF').gt.0) then
4353 else if (index(controlcard,'MAP').gt.0) then
4355 call readi(controlcard,'MAP',nmap,0)
4356 else if (index(controlcard,'CSA').gt.0) then
4358 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4360 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4363 !fcm else if (index(controlcard,'MCMF').gt.0) then
4365 else if (index(controlcard,'SOFTREG').gt.0) then
4367 else if (index(controlcard,'CHECK_BOND').gt.0) then
4369 else if (index(controlcard,'TEST').gt.0) then
4371 else if (index(controlcard,'MD').gt.0) then
4373 else if (index(controlcard,'RE ').gt.0) then
4377 lmuca=index(controlcard,'MUCA').gt.0
4378 call readi(controlcard,'MUCADYN',mucadyn,0)
4379 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4380 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4382 write (iout,*) 'MUCADYN=',mucadyn
4383 write (iout,*) 'MUCASMOOTH=',muca_smooth
4386 iscode=index(controlcard,'ONE_LETTER')
4387 indphi=index(controlcard,'PHI')
4388 indback=index(controlcard,'BACK')
4389 iranconf=index(controlcard,'RAND_CONF')
4390 i2ndstr=index(controlcard,'USE_SEC_PRED')
4391 gradout=index(controlcard,'GRADOUT').gt.0
4392 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4393 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4394 if (me.eq.king .or. .not.out1file ) &
4395 write (iout,*) "DISTCHAINMAX",distchainmax
4397 if(me.eq.king.or..not.out1file) &
4398 write (iout,'(2a)') diagmeth(kdiag),&
4399 ' routine used to diagonalize matrices.'
4400 if (shield_mode.gt.0) then
4401 pi=4.0D0*datan(1.0D0)
4402 !C VSolvSphere the volume of solving sphere
4404 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4405 !C there will be no distinction between proline peptide group and normal peptide
4406 !C group in case of shielding parameters
4407 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4408 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4409 write (iout,*) VSolvSphere,VSolvSphere_div
4410 !C long axis of side chain
4412 ! long_r_sidechain(i)=vbldsc0(1,i)
4413 ! short_r_sidechain(i)=sigma0(i)
4414 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4420 end subroutine read_control
4421 !-----------------------------------------------------------------------------
4422 subroutine read_REMDpar
4424 ! Read REMD settings
4431 use control_data, only:out1file
4432 ! implicit real*8 (a-h,o-z)
4433 ! include 'DIMENSIONS'
4434 ! include 'COMMON.IOUNITS'
4435 ! include 'COMMON.TIME1'
4436 ! include 'COMMON.MD'
4439 !el include 'COMMON.LANGEVIN'
4441 !el include 'COMMON.LANGEVIN.lang0'
4443 ! include 'COMMON.INTERACT'
4444 ! include 'COMMON.NAMES'
4445 ! include 'COMMON.GEO'
4446 ! include 'COMMON.REMD'
4447 ! include 'COMMON.CONTROL'
4448 ! include 'COMMON.SETUP'
4449 ! character(len=80) :: ucase
4450 character(len=320) :: controlcard
4451 character(len=3200) :: controlcard1
4452 integer :: iremd_m_total
4455 ! real(kind=8) :: var,ene
4457 if(me.eq.king.or..not.out1file) &
4458 write (iout,*) "REMD setup"
4460 call card_concat(controlcard,.true.)
4461 call readi(controlcard,"NREP",nrep,3)
4462 call readi(controlcard,"NSTEX",nstex,1000)
4463 call reada(controlcard,"RETMIN",retmin,10.0d0)
4464 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4465 mremdsync=(index(controlcard,'SYNC').gt.0)
4466 call readi(controlcard,"NSYN",i_sync_step,100)
4467 restart1file=(index(controlcard,'REST1FILE').gt.0)
4468 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4469 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4470 if(max_cache_traj_use.gt.max_cache_traj) &
4471 max_cache_traj_use=max_cache_traj
4472 if(me.eq.king.or..not.out1file) then
4473 !d if (traj1file) then
4474 !rc caching is in testing - NTWX is not ignored
4475 !d write (iout,*) "NTWX value is ignored"
4476 !d write (iout,*) " trajectory is stored to one file by master"
4477 !d write (iout,*) " before exchange at NSTEX intervals"
4479 write (iout,*) "NREP= ",nrep
4480 write (iout,*) "NSTEX= ",nstex
4481 write (iout,*) "SYNC= ",mremdsync
4482 write (iout,*) "NSYN= ",i_sync_step
4483 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4486 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4487 if (index(controlcard,'TLIST').gt.0) then
4489 call card_concat(controlcard1,.true.)
4490 read(controlcard1,*) (remd_t(i),i=1,nrep)
4491 if(me.eq.king.or..not.out1file) &
4492 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4495 if (index(controlcard,'MLIST').gt.0) then
4497 call card_concat(controlcard1,.true.)
4498 read(controlcard1,*) (remd_m(i),i=1,nrep)
4499 if(me.eq.king.or..not.out1file) then
4500 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4503 iremd_m_total=iremd_m_total+remd_m(i)
4505 write (iout,*) 'Total number of replicas ',iremd_m_total
4508 if(me.eq.king.or..not.out1file) &
4509 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4511 end subroutine read_REMDpar
4512 !-----------------------------------------------------------------------------
4513 subroutine read_MDpar
4517 use control_data, only: r_cut,rlamb,out1file
4519 use geometry_data, only: pi
4521 ! implicit real*8 (a-h,o-z)
4522 ! include 'DIMENSIONS'
4523 ! include 'COMMON.IOUNITS'
4524 ! include 'COMMON.TIME1'
4525 ! include 'COMMON.MD'
4528 !el include 'COMMON.LANGEVIN'
4530 !el include 'COMMON.LANGEVIN.lang0'
4532 ! include 'COMMON.INTERACT'
4533 ! include 'COMMON.NAMES'
4534 ! include 'COMMON.GEO'
4535 ! include 'COMMON.SETUP'
4536 ! include 'COMMON.CONTROL'
4537 ! include 'COMMON.SPLITELE'
4538 ! character(len=80) :: ucase
4539 character(len=320) :: controlcard
4544 call card_concat(controlcard,.true.)
4545 call readi(controlcard,"NSTEP",n_timestep,1000000)
4546 call readi(controlcard,"NTWE",ntwe,100)
4547 call readi(controlcard,"NTWX",ntwx,1000)
4548 call reada(controlcard,"DT",d_time,1.0d-1)
4549 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4550 call reada(controlcard,"DAMAX",damax,1.0d1)
4551 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4552 call readi(controlcard,"LANG",lang,0)
4553 RESPA = index(controlcard,"RESPA") .gt. 0
4554 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4555 ntime_split0=ntime_split
4556 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4557 ntime_split0=ntime_split
4558 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4559 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4560 rest = index(controlcard,"REST").gt.0
4561 tbf = index(controlcard,"TBF").gt.0
4562 usampl = index(controlcard,"USAMPL").gt.0
4563 mdpdb = index(controlcard,"MDPDB").gt.0
4564 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4565 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4566 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4567 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4568 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4569 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4570 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4571 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4572 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4573 large = index(controlcard,"LARGE").gt.0
4574 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4575 rattle = index(controlcard,"RATTLE").gt.0
4576 preminim=(index(controlcard,'PREMINIM').gt.0)
4577 write (iout,*) "PREMINIM ",preminim
4578 dccart=(index(controlcard,'CART').gt.0)
4579 if (preminim) call read_minim
4580 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4586 if(me.eq.king.or..not.out1file) then
4588 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4590 write (iout,'(a)') "The units are:"
4591 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4592 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4593 " acceleration: angstrom/(48.9 fs)**2"
4594 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4596 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4597 write (iout,'(a60,f10.5,a)') &
4598 "Initial time step of numerical integration:",d_time,&
4600 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4602 write (iout,'(2a,i4,a)') &
4603 "A-MTS algorithm used; initial time step for fast-varying",&
4604 " short-range forces split into",ntime_split," steps."
4605 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4606 r_cut," lambda",rlamb
4608 write (iout,'(2a,f10.5)') &
4609 "Maximum acceleration threshold to reduce the time step",&
4610 "/increase split number:",damax
4611 write (iout,'(2a,f10.5)') &
4612 "Maximum predicted energy drift to reduce the timestep",&
4613 "/increase split number:",edriftmax
4614 write (iout,'(a60,f10.5)') &
4615 "Maximum velocity threshold to reduce velocities:",dvmax
4616 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4617 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4618 if (rattle) write (iout,'(a60)') &
4619 "Rattle algorithm used to constrain the virtual bonds"
4623 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4624 call reada(controlcard,"RWAT",rwat,1.4d0)
4625 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4626 surfarea=index(controlcard,"SURFAREA").gt.0
4627 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4628 if(me.eq.king.or..not.out1file)then
4629 write (iout,'(/a,$)') "Langevin dynamics calculation"
4631 write (iout,'(a/)') &
4632 " with direct integration of Langevin equations"
4633 else if (lang.eq.2) then
4634 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4635 else if (lang.eq.3) then
4636 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4637 else if (lang.eq.4) then
4638 write (iout,'(a/)') " in overdamped mode"
4640 write (iout,'(//a,i5)') &
4641 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4644 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4645 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4646 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4647 write (iout,'(a60,f10.5)') &
4648 "Scaling factor of the friction forces:",scal_fric
4649 if (surfarea) write (iout,'(2a,i10,a)') &
4650 "Friction coefficients will be scaled by solvent-accessible",&
4651 " surface area every",reset_fricmat," steps."
4653 ! Calculate friction coefficients and bounds of stochastic forces
4654 eta=6*pi*cPoise*etawat
4655 if(me.eq.king.or..not.out1file) &
4656 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4659 do j=1,5 !types of molecules
4660 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4661 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4663 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4664 do j=1,5 !types of molecules
4666 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4667 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4671 if(me.eq.king.or..not.out1file)then
4672 write (iout,'(/2a/)') &
4673 "Radii of site types and friction coefficients and std's of",&
4674 " stochastic forces of fully exposed sites"
4675 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4677 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4678 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4682 if(me.eq.king.or..not.out1file)then
4683 write (iout,'(a)') "Berendsen bath calculation"
4684 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4685 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4687 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4688 count_reset_moment," steps"
4690 write (iout,'(a,i10,a)') &
4691 "Velocities will be reset at random every",count_reset_vel,&
4695 if(me.eq.king.or..not.out1file) &
4696 write (iout,'(a31)') "Microcanonical mode calculation"
4698 if(me.eq.king.or..not.out1file)then
4699 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4701 write(iout,*) "MD running with constraints."
4702 write(iout,*) "Equilibration time ", eq_time, " mtus."
4703 write(iout,*) "Constraining ", nfrag," fragments."
4704 write(iout,*) "Length of each fragment, weight and q0:"
4706 write (iout,*) "Set of restraints #",iset
4708 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4709 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4711 write(iout,*) "constraints between ", npair, "fragments."
4712 write(iout,*) "constraint pairs, weights and q0:"
4714 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4715 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4717 write(iout,*) "angle constraints within ", nfrag_back,&
4718 "backbone fragments."
4719 write(iout,*) "fragment, weights:"
4721 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4722 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4723 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4726 iset=mod(kolor,nset)+1
4729 if(me.eq.king.or..not.out1file) &
4730 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4732 end subroutine read_MDpar
4733 !-----------------------------------------------------------------------------
4737 ! implicit real*8 (a-h,o-z)
4738 ! include 'DIMENSIONS'
4739 ! include 'COMMON.MAP'
4740 ! include 'COMMON.IOUNITS'
4741 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4742 character(len=80) :: mapcard !,ucase
4745 ! real(kind=8) :: var,ene
4748 read (inp,'(a)') mapcard
4749 mapcard=ucase(mapcard)
4750 if (index(mapcard,'PHI').gt.0) then
4752 else if (index(mapcard,'THE').gt.0) then
4754 else if (index(mapcard,'ALP').gt.0) then
4756 else if (index(mapcard,'OME').gt.0) then
4759 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4760 stop 'Error - illegal variable spec in MAP card.'
4762 call readi (mapcard,'RES1',res1(imap),0)
4763 call readi (mapcard,'RES2',res2(imap),0)
4764 if (res1(imap).eq.0) then
4765 res1(imap)=res2(imap)
4766 else if (res2(imap).eq.0) then
4767 res2(imap)=res1(imap)
4769 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4770 write (iout,'(a)') &
4771 'Error - illegal definition of variable group in MAP.'
4772 stop 'Error - illegal definition of variable group in MAP.'
4774 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4775 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4776 call readi(mapcard,'NSTEP',nstep(imap),0)
4777 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4778 write (iout,'(a)') &
4779 'Illegal boundary and/or step size specification in MAP.'
4780 stop 'Illegal boundary and/or step size specification in MAP.'
4784 end subroutine map_read
4785 !-----------------------------------------------------------------------------
4788 use control_data, only: vdisulf
4790 ! implicit real*8 (a-h,o-z)
4791 ! include 'DIMENSIONS'
4792 ! include 'COMMON.IOUNITS'
4793 ! include 'COMMON.GEO'
4794 ! include 'COMMON.CSA'
4795 ! include 'COMMON.BANK'
4796 ! include 'COMMON.CONTROL'
4797 ! character(len=80) :: ucase
4798 character(len=620) :: mcmcard
4800 ! integer :: ntf,ik,iw_pdb
4801 ! real(kind=8) :: var,ene
4803 call card_concat(mcmcard,.true.)
4805 call readi(mcmcard,'NCONF',nconf,50)
4806 call readi(mcmcard,'NADD',nadd,0)
4807 call readi(mcmcard,'JSTART',jstart,1)
4808 call readi(mcmcard,'JEND',jend,1)
4809 call readi(mcmcard,'NSTMAX',nstmax,500000)
4810 call readi(mcmcard,'N0',n0,1)
4811 call readi(mcmcard,'N1',n1,6)
4812 call readi(mcmcard,'N2',n2,4)
4813 call readi(mcmcard,'N3',n3,0)
4814 call readi(mcmcard,'N4',n4,0)
4815 call readi(mcmcard,'N5',n5,0)
4816 call readi(mcmcard,'N6',n6,10)
4817 call readi(mcmcard,'N7',n7,0)
4818 call readi(mcmcard,'N8',n8,0)
4819 call readi(mcmcard,'N9',n9,0)
4820 call readi(mcmcard,'N14',n14,0)
4821 call readi(mcmcard,'N15',n15,0)
4822 call readi(mcmcard,'N16',n16,0)
4823 call readi(mcmcard,'N17',n17,0)
4824 call readi(mcmcard,'N18',n18,0)
4826 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4828 call readi(mcmcard,'NDIFF',ndiff,2)
4829 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4830 call readi(mcmcard,'IS1',is1,1)
4831 call readi(mcmcard,'IS2',is2,8)
4832 call readi(mcmcard,'NRAN0',nran0,4)
4833 call readi(mcmcard,'NRAN1',nran1,2)
4834 call readi(mcmcard,'IRR',irr,1)
4835 call readi(mcmcard,'NSEED',nseed,20)
4836 call readi(mcmcard,'NTOTAL',ntotal,10000)
4837 call reada(mcmcard,'CUT1',cut1,2.0d0)
4838 call reada(mcmcard,'CUT2',cut2,5.0d0)
4839 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4840 call readi(mcmcard,'ICMAX',icmax,3)
4841 call readi(mcmcard,'IRESTART',irestart,0)
4842 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4845 call reada(mcmcard,'DELE',dele,20.0d0)
4846 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4847 call readi(mcmcard,'IREF',iref,0)
4848 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4849 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4850 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4851 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4852 write (iout,*) "NCONF_IN",nconf_in
4854 end subroutine csaread
4855 !-----------------------------------------------------------------------------
4859 use control_data, only: MaxMoveType
4862 ! implicit real*8 (a-h,o-z)
4863 ! include 'DIMENSIONS'
4864 ! include 'COMMON.MCM'
4865 ! include 'COMMON.MCE'
4866 ! include 'COMMON.IOUNITS'
4867 ! character(len=80) :: ucase
4868 character(len=320) :: mcmcard
4871 ! real(kind=8) :: var,ene
4873 call card_concat(mcmcard,.true.)
4874 call readi(mcmcard,'MAXACC',maxacc,100)
4875 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4876 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4877 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4878 call readi(mcmcard,'MAXREPM',maxrepm,200)
4879 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4880 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4881 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4882 call reada(mcmcard,'E_UP',e_up,5.0D0)
4883 call reada(mcmcard,'DELTE',delte,0.1D0)
4884 call readi(mcmcard,'NSWEEP',nsweep,5)
4885 call readi(mcmcard,'NSTEPH',nsteph,0)
4886 call readi(mcmcard,'NSTEPC',nstepc,0)
4887 call reada(mcmcard,'TMIN',tmin,298.0D0)
4888 call reada(mcmcard,'TMAX',tmax,298.0D0)
4889 call readi(mcmcard,'NWINDOW',nwindow,0)
4890 call readi(mcmcard,'PRINT_MC',print_mc,0)
4891 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4892 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4893 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4894 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4895 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4896 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4897 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4898 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4899 if (nwindow.gt.0) then
4900 allocate(winstart(nwindow)) !!el (maxres)
4901 allocate(winend(nwindow)) !!el
4902 allocate(winlen(nwindow)) !!el
4903 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4905 winlen(i)=winend(i)-winstart(i)+1
4908 if (tmax.lt.tmin) tmax=tmin
4909 if (tmax.eq.tmin) then
4913 if (nstepc.gt.0 .and. nsteph.gt.0) then
4914 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4915 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4917 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4918 ! Probabilities of different move types
4919 sumpro_type(0)=0.0D0
4920 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4921 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4922 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4923 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4924 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4925 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4926 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4928 print *,'i',i,' sumprotype',sumpro_type(i)
4929 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4930 print *,'i',i,' sumprotype',sumpro_type(i)
4933 end subroutine mcmread
4934 !-----------------------------------------------------------------------------
4935 subroutine read_minim
4938 ! implicit real*8 (a-h,o-z)
4939 ! include 'DIMENSIONS'
4940 ! include 'COMMON.MINIM'
4941 ! include 'COMMON.IOUNITS'
4942 ! character(len=80) :: ucase
4943 character(len=320) :: minimcard
4945 ! integer :: ntf,ik,iw_pdb
4946 ! real(kind=8) :: var,ene
4948 call card_concat(minimcard,.true.)
4949 call readi(minimcard,'MAXMIN',maxmin,2000)
4950 call readi(minimcard,'MAXFUN',maxfun,5000)
4951 call readi(minimcard,'MINMIN',minmin,maxmin)
4952 call readi(minimcard,'MINFUN',minfun,maxmin)
4953 call reada(minimcard,'TOLF',tolf,1.0D-2)
4954 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4955 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4956 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4957 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4958 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4959 'Options in energy minimization:'
4960 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4961 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4962 'MinMin:',MinMin,' MinFun:',MinFun,&
4963 ' TolF:',TolF,' RTolF:',RTolF
4965 end subroutine read_minim
4966 !-----------------------------------------------------------------------------
4967 subroutine openunits
4969 use MD_data, only: usampl
4972 use control_data, only:out1file
4973 use control, only: getenv_loc
4974 ! implicit real*8 (a-h,o-z)
4975 ! include 'DIMENSIONS'
4978 character(len=16) :: form,nodename
4979 integer :: nodelen,ierror,npos
4981 ! include 'COMMON.SETUP'
4982 ! include 'COMMON.IOUNITS'
4983 ! include 'COMMON.MD'
4984 ! include 'COMMON.CONTROL'
4985 integer :: lenpre,lenpot,lentmp !,ilen
4987 character(len=3) :: out1file_text !,ucase
4988 character(len=3) :: ll
4991 ! integer :: ntf,ik,iw_pdb
4992 ! real(kind=8) :: var,ene
4994 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4995 call getenv_loc("PREFIX",prefix)
4997 call getenv_loc("POT",pot)
4998 call getenv_loc("DIRTMP",tmpdir)
4999 call getenv_loc("CURDIR",curdir)
5000 call getenv_loc("OUT1FILE",out1file_text)
5001 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5002 out1file_text=ucase(out1file_text)
5003 if (out1file_text(1:1).eq."Y") then
5006 out1file=fg_rank.gt.0
5011 if (lentmp.gt.0) then
5012 write (*,'(80(1h!))')
5013 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5014 write (*,'(80(1h!))')
5015 write (*,*)"All output files will be on node /tmp directory."
5017 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5018 if (me.eq.king) then
5019 write (*,*) "The master node is ",nodename
5020 else if (fg_rank.eq.0) then
5021 write (*,*) "I am the CG slave node ",nodename
5023 write (*,*) "I am the FG slave node ",nodename
5026 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5027 lenpre = lentmp+lenpre+1
5029 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5030 ! Get the names and open the input files
5031 #if defined(WINIFL) || defined(WINPGI)
5032 open(1,file=pref_orig(:ilen(pref_orig))// &
5033 '.inp',status='old',readonly,shared)
5034 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5035 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5036 ! Get parameter filenames and open the parameter files.
5037 call getenv_loc('BONDPAR',bondname)
5038 open (ibond,file=bondname,status='old',readonly,shared)
5039 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5040 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5041 call getenv_loc('THETPAR',thetname)
5042 open (ithep,file=thetname,status='old',readonly,shared)
5043 call getenv_loc('ROTPAR',rotname)
5044 open (irotam,file=rotname,status='old',readonly,shared)
5045 call getenv_loc('TORPAR',torname)
5046 open (itorp,file=torname,status='old',readonly,shared)
5047 call getenv_loc('TORDPAR',tordname)
5048 open (itordp,file=tordname,status='old',readonly,shared)
5049 call getenv_loc('FOURIER',fouriername)
5050 open (ifourier,file=fouriername,status='old',readonly,shared)
5051 call getenv_loc('ELEPAR',elename)
5052 open (ielep,file=elename,status='old',readonly,shared)
5053 call getenv_loc('SIDEPAR',sidename)
5054 open (isidep,file=sidename,status='old',readonly,shared)
5056 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5057 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5058 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5059 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5060 call getenv_loc('TORPAR_NUCL',torname_nucl)
5061 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5062 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5063 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5064 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5065 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5068 #elif (defined CRAY) || (defined AIX)
5069 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5071 ! print *,"Processor",myrank," opened file 1"
5072 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5073 ! print *,"Processor",myrank," opened file 9"
5074 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5075 ! Get parameter filenames and open the parameter files.
5076 call getenv_loc('BONDPAR',bondname)
5077 open (ibond,file=bondname,status='old',action='read')
5078 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5079 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5081 ! print *,"Processor",myrank," opened file IBOND"
5082 call getenv_loc('THETPAR',thetname)
5083 open (ithep,file=thetname,status='old',action='read')
5084 ! print *,"Processor",myrank," opened file ITHEP"
5085 call getenv_loc('ROTPAR',rotname)
5086 open (irotam,file=rotname,status='old',action='read')
5087 ! print *,"Processor",myrank," opened file IROTAM"
5088 call getenv_loc('TORPAR',torname)
5089 open (itorp,file=torname,status='old',action='read')
5090 ! print *,"Processor",myrank," opened file ITORP"
5091 call getenv_loc('TORDPAR',tordname)
5092 open (itordp,file=tordname,status='old',action='read')
5093 ! print *,"Processor",myrank," opened file ITORDP"
5094 call getenv_loc('SCCORPAR',sccorname)
5095 open (isccor,file=sccorname,status='old',action='read')
5096 ! print *,"Processor",myrank," opened file ISCCOR"
5097 call getenv_loc('FOURIER',fouriername)
5098 open (ifourier,file=fouriername,status='old',action='read')
5099 ! print *,"Processor",myrank," opened file IFOURIER"
5100 call getenv_loc('ELEPAR',elename)
5101 open (ielep,file=elename,status='old',action='read')
5102 ! print *,"Processor",myrank," opened file IELEP"
5103 call getenv_loc('SIDEPAR',sidename)
5104 open (isidep,file=sidename,status='old',action='read')
5106 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5107 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5108 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5109 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5110 call getenv_loc('TORPAR_NUCL',torname_nucl)
5111 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5112 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5113 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5114 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5115 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5117 call getenv_loc('LIPTRANPAR',liptranname)
5118 open (iliptranpar,file=liptranname,status='old',action='read')
5119 call getenv_loc('TUBEPAR',tubename)
5120 open (itube,file=tubename,status='old',action='read')
5121 call getenv_loc('IONPAR',ionname)
5122 open (iion,file=ionname,status='old',action='read')
5124 ! print *,"Processor",myrank," opened file ISIDEP"
5125 ! print *,"Processor",myrank," opened parameter files"
5127 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5128 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5129 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5130 ! Get parameter filenames and open the parameter files.
5131 call getenv_loc('BONDPAR',bondname)
5132 open (ibond,file=bondname,status='old')
5133 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5134 open (ibond_nucl,file=bondname_nucl,status='old')
5136 call getenv_loc('THETPAR',thetname)
5137 open (ithep,file=thetname,status='old')
5138 call getenv_loc('ROTPAR',rotname)
5139 open (irotam,file=rotname,status='old')
5140 call getenv_loc('TORPAR',torname)
5141 open (itorp,file=torname,status='old')
5142 call getenv_loc('TORDPAR',tordname)
5143 open (itordp,file=tordname,status='old')
5144 call getenv_loc('SCCORPAR',sccorname)
5145 open (isccor,file=sccorname,status='old')
5146 call getenv_loc('FOURIER',fouriername)
5147 open (ifourier,file=fouriername,status='old')
5148 call getenv_loc('ELEPAR',elename)
5149 open (ielep,file=elename,status='old')
5150 call getenv_loc('SIDEPAR',sidename)
5151 open (isidep,file=sidename,status='old')
5153 open (ithep_nucl,file=thetname_nucl,status='old')
5154 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5155 open (irotam_nucl,file=rotname_nucl,status='old')
5156 call getenv_loc('TORPAR_NUCL',torname_nucl)
5157 open (itorp_nucl,file=torname_nucl,status='old')
5158 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5159 open (itordp_nucl,file=tordname_nucl,status='old')
5160 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5161 open (isidep_nucl,file=sidename_nucl,status='old')
5163 call getenv_loc('LIPTRANPAR',liptranname)
5164 open (iliptranpar,file=liptranname,status='old')
5165 call getenv_loc('TUBEPAR',tubename)
5166 open (itube,file=tubename,status='old')
5167 call getenv_loc('IONPAR',ionname)
5168 open (iion,file=ionname,status='old')
5170 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5172 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5173 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5174 ! Get parameter filenames and open the parameter files.
5175 call getenv_loc('BONDPAR',bondname)
5176 open (ibond,file=bondname,status='old',action='read')
5177 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5178 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5179 call getenv_loc('THETPAR',thetname)
5180 open (ithep,file=thetname,status='old',action='read')
5181 call getenv_loc('ROTPAR',rotname)
5182 open (irotam,file=rotname,status='old',action='read')
5183 call getenv_loc('TORPAR',torname)
5184 open (itorp,file=torname,status='old',action='read')
5185 call getenv_loc('TORDPAR',tordname)
5186 open (itordp,file=tordname,status='old',action='read')
5187 call getenv_loc('SCCORPAR',sccorname)
5188 open (isccor,file=sccorname,status='old',action='read')
5190 call getenv_loc('THETPARPDB',thetname_pdb)
5191 print *,"thetname_pdb ",thetname_pdb
5192 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5193 print *,ithep_pdb," opened"
5195 call getenv_loc('FOURIER',fouriername)
5196 open (ifourier,file=fouriername,status='old',readonly)
5197 call getenv_loc('ELEPAR',elename)
5198 open (ielep,file=elename,status='old',readonly)
5199 call getenv_loc('SIDEPAR',sidename)
5200 open (isidep,file=sidename,status='old',readonly)
5202 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5203 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5204 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5205 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5206 call getenv_loc('TORPAR_NUCL',torname_nucl)
5207 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5208 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5209 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5210 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5211 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5212 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5213 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5214 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5215 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5216 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5217 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5218 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5219 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5222 call getenv_loc('LIPTRANPAR',liptranname)
5223 open (iliptranpar,file=liptranname,status='old',action='read')
5224 call getenv_loc('TUBEPAR',tubename)
5225 open (itube,file=tubename,status='old',action='read')
5226 call getenv_loc('IONPAR',ionname)
5227 open (iion,file=ionname,status='old',action='read')
5230 call getenv_loc('ROTPARPDB',rotname_pdb)
5231 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5234 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5235 #if defined(WINIFL) || defined(WINPGI)
5236 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5237 #elif (defined CRAY) || (defined AIX)
5238 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5240 open (iscpp_nucl,file=scpname_nucl,status='old')
5242 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5247 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5248 ! Use -DOLDSCP to use hard-coded constants instead.
5250 call getenv_loc('SCPPAR',scpname)
5251 #if defined(WINIFL) || defined(WINPGI)
5252 open (iscpp,file=scpname,status='old',readonly,shared)
5253 #elif (defined CRAY) || (defined AIX)
5254 open (iscpp,file=scpname,status='old',action='read')
5256 open (iscpp,file=scpname,status='old')
5258 open (iscpp,file=scpname,status='old',action='read')
5261 call getenv_loc('PATTERN',patname)
5262 #if defined(WINIFL) || defined(WINPGI)
5263 open (icbase,file=patname,status='old',readonly,shared)
5264 #elif (defined CRAY) || (defined AIX)
5265 open (icbase,file=patname,status='old',action='read')
5267 open (icbase,file=patname,status='old')
5269 open (icbase,file=patname,status='old',action='read')
5272 ! Open output file only for CG processes
5273 ! print *,"Processor",myrank," fg_rank",fg_rank
5274 if (fg_rank.eq.0) then
5276 if (nodes.eq.1) then
5279 npos = dlog10(dfloat(nodes-1))+1
5281 if (npos.lt.3) npos=3
5282 write (liczba,'(i1)') npos
5283 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5285 write (liczba,form) me
5286 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5287 liczba(:ilen(liczba))
5288 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5290 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5292 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5293 liczba(:ilen(liczba))//'.mol2'
5294 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5295 liczba(:ilen(liczba))//'.stat'
5297 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5298 //liczba(:ilen(liczba))//'.stat')
5299 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5302 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5303 liczba(:ilen(liczba))//'.const'
5308 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5309 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5310 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5311 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5312 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5314 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5316 rest2name=prefix(:ilen(prefix))//'.rst'
5318 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5321 #if defined(AIX) || defined(PGI)
5322 if (me.eq.king .or. .not. out1file) &
5323 open(iout,file=outname,status='unknown')
5325 if (fg_rank.gt.0) then
5326 write (liczba,'(i3.3)') myrank/nfgtasks
5327 write (ll,'(bz,i3.3)') fg_rank
5328 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5333 open(igeom,file=intname,status='unknown',position='append')
5334 open(ipdb,file=pdbname,status='unknown')
5335 open(imol2,file=mol2name,status='unknown')
5336 open(istat,file=statname,status='unknown',position='append')
5338 !1out open(iout,file=outname,status='unknown')
5341 if (me.eq.king .or. .not.out1file) &
5342 open(iout,file=outname,status='unknown')
5344 if (fg_rank.gt.0) then
5345 write (liczba,'(i3.3)') myrank/nfgtasks
5346 write (ll,'(bz,i3.3)') fg_rank
5347 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5352 open(igeom,file=intname,status='unknown',access='append')
5353 open(ipdb,file=pdbname,status='unknown')
5354 open(imol2,file=mol2name,status='unknown')
5355 open(istat,file=statname,status='unknown',access='append')
5357 !1out open(iout,file=outname,status='unknown')
5360 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5361 csa_seed=prefix(:lenpre)//'.CSA.seed'
5362 csa_history=prefix(:lenpre)//'.CSA.history'
5363 csa_bank=prefix(:lenpre)//'.CSA.bank'
5364 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5365 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5366 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5367 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5368 csa_int=prefix(:lenpre)//'.int'
5369 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5370 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5371 csa_in=prefix(:lenpre)//'.CSA.in'
5372 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5375 write (iout,'(80(1h-))')
5376 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5377 write (iout,'(80(1h-))')
5378 write (iout,*) "Input file : ",&
5379 pref_orig(:ilen(pref_orig))//'.inp'
5380 write (iout,*) "Output file : ",&
5381 outname(:ilen(outname))
5383 write (iout,*) "Sidechain potential file : ",&
5384 sidename(:ilen(sidename))
5386 write (iout,*) "SCp potential file : ",&
5387 scpname(:ilen(scpname))
5389 write (iout,*) "Electrostatic potential file : ",&
5390 elename(:ilen(elename))
5391 write (iout,*) "Cumulant coefficient file : ",&
5392 fouriername(:ilen(fouriername))
5393 write (iout,*) "Torsional parameter file : ",&
5394 torname(:ilen(torname))
5395 write (iout,*) "Double torsional parameter file : ",&
5396 tordname(:ilen(tordname))
5397 write (iout,*) "SCCOR parameter file : ",&
5398 sccorname(:ilen(sccorname))
5399 write (iout,*) "Bond & inertia constant file : ",&
5400 bondname(:ilen(bondname))
5401 write (iout,*) "Bending parameter file : ",&
5402 thetname(:ilen(thetname))
5403 write (iout,*) "Rotamer parameter file : ",&
5404 rotname(:ilen(rotname))
5407 write (iout,*) "Thetpdb parameter file : ",&
5408 thetname_pdb(:ilen(thetname_pdb))
5411 write (iout,*) "Threading database : ",&
5412 patname(:ilen(patname))
5414 write (iout,*)" DIRTMP : ",&
5416 write (iout,'(80(1h-))')
5419 end subroutine openunits
5420 !-----------------------------------------------------------------------------
5423 use geometry_data, only: nres,dc
5425 ! implicit real*8 (a-h,o-z)
5426 ! include 'DIMENSIONS'
5427 ! include 'COMMON.CHAIN'
5428 ! include 'COMMON.IOUNITS'
5429 ! include 'COMMON.MD'
5432 ! real(kind=8) :: var,ene
5434 open(irest2,file=rest2name,status='unknown')
5435 read(irest2,*) totT,EK,potE,totE,t_bath
5438 ! AL 4/17/17: Now reading d_t(0,:) too
5440 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5443 ! AL 4/17/17: Now reading d_c(0,:) too
5445 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5448 read (irest2,*) iset
5452 end subroutine readrst
5453 !-----------------------------------------------------------------------------
5454 subroutine read_fragments
5458 use control_data, only:out1file
5461 ! implicit real*8 (a-h,o-z)
5462 ! include 'DIMENSIONS'
5466 ! include 'COMMON.SETUP'
5467 ! include 'COMMON.CHAIN'
5468 ! include 'COMMON.IOUNITS'
5469 ! include 'COMMON.MD'
5470 ! include 'COMMON.CONTROL'
5473 ! real(kind=8) :: var,ene
5475 read(inp,*) nset,nfrag,npair,nfrag_back
5477 !el from module energy
5478 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5479 if(.not.allocated(wfrag_back)) then
5480 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5481 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5483 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5484 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5486 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5487 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5490 if(me.eq.king.or..not.out1file) &
5491 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5492 " nfrag_back",nfrag_back
5494 read(inp,*) mset(iset)
5496 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5498 if(me.eq.king.or..not.out1file) &
5499 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5500 ifrag(2,i,iset), qinfrag(i,iset)
5503 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5505 if(me.eq.king.or..not.out1file) &
5506 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5507 ipair(2,i,iset), qinpair(i,iset)
5510 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5511 wfrag_back(3,i,iset),&
5512 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5513 if(me.eq.king.or..not.out1file) &
5514 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5515 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5519 end subroutine read_fragments
5520 !-----------------------------------------------------------------------------
5522 !-----------------------------------------------------------------------------
5526 ! implicit real*8 (a-h,o-z)
5527 ! include 'DIMENSIONS'
5528 ! include 'COMMON.CSA'
5529 ! include 'COMMON.BANK'
5530 ! include 'COMMON.IOUNITS'
5532 ! integer :: ntf,ik,iw_pdb
5533 ! real(kind=8) :: var,ene
5535 open(icsa_in,file=csa_in,status="old",err=100)
5536 read(icsa_in,*) nconf
5537 read(icsa_in,*) jstart,jend
5538 read(icsa_in,*) nstmax
5539 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5540 read(icsa_in,*) nran0,nran1,irr
5541 read(icsa_in,*) nseed
5542 read(icsa_in,*) ntotal,cut1,cut2
5543 read(icsa_in,*) estop
5544 read(icsa_in,*) icmax,irestart
5545 read(icsa_in,*) ntbankm,dele,difcut
5546 read(icsa_in,*) iref,rmscut,pnccut
5547 read(icsa_in,*) ndiff
5554 end subroutine csa_read
5555 !-----------------------------------------------------------------------------
5556 subroutine initial_write
5559 ! implicit real*8 (a-h,o-z)
5560 ! include 'DIMENSIONS'
5561 ! include 'COMMON.CSA'
5562 ! include 'COMMON.BANK'
5563 ! include 'COMMON.IOUNITS'
5565 ! integer :: ntf,ik,iw_pdb
5566 ! real(kind=8) :: var,ene
5568 open(icsa_seed,file=csa_seed,status="unknown")
5569 write(icsa_seed,*) "seed"
5571 #if defined(AIX) || defined(PGI)
5572 open(icsa_history,file=csa_history,status="unknown",&
5575 open(icsa_history,file=csa_history,status="unknown",&
5578 write(icsa_history,*) nconf
5579 write(icsa_history,*) jstart,jend
5580 write(icsa_history,*) nstmax
5581 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5582 write(icsa_history,*) nran0,nran1,irr
5583 write(icsa_history,*) nseed
5584 write(icsa_history,*) ntotal,cut1,cut2
5585 write(icsa_history,*) estop
5586 write(icsa_history,*) icmax,irestart
5587 write(icsa_history,*) ntbankm,dele,difcut
5588 write(icsa_history,*) iref,rmscut,pnccut
5589 write(icsa_history,*) ndiff
5591 write(icsa_history,*)
5594 open(icsa_bank1,file=csa_bank1,status="unknown")
5595 write(icsa_bank1,*) 0
5599 end subroutine initial_write
5600 !-----------------------------------------------------------------------------
5601 subroutine restart_write
5604 ! implicit real*8 (a-h,o-z)
5605 ! include 'DIMENSIONS'
5606 ! include 'COMMON.IOUNITS'
5607 ! include 'COMMON.CSA'
5608 ! include 'COMMON.BANK'
5610 ! integer :: ntf,ik,iw_pdb
5611 ! real(kind=8) :: var,ene
5613 #if defined(AIX) || defined(PGI)
5614 open(icsa_history,file=csa_history,position="append")
5616 open(icsa_history,file=csa_history,access="append")
5618 write(icsa_history,*)
5619 write(icsa_history,*) "This is restart"
5620 write(icsa_history,*)
5621 write(icsa_history,*) nconf
5622 write(icsa_history,*) jstart,jend
5623 write(icsa_history,*) nstmax
5624 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5625 write(icsa_history,*) nran0,nran1,irr
5626 write(icsa_history,*) nseed
5627 write(icsa_history,*) ntotal,cut1,cut2
5628 write(icsa_history,*) estop
5629 write(icsa_history,*) icmax,irestart
5630 write(icsa_history,*) ntbankm,dele,difcut
5631 write(icsa_history,*) iref,rmscut,pnccut
5632 write(icsa_history,*) ndiff
5633 write(icsa_history,*)
5634 write(icsa_history,*) "irestart is: ", irestart
5636 write(icsa_history,*)
5640 end subroutine restart_write
5641 !-----------------------------------------------------------------------------
5643 !-----------------------------------------------------------------------------
5644 subroutine write_pdb(npdb,titelloc,ee)
5646 ! implicit real*8 (a-h,o-z)
5647 ! include 'DIMENSIONS'
5648 ! include 'COMMON.IOUNITS'
5649 character(len=50) :: titelloc1
5650 character*(*) :: titelloc
5651 character(len=3) :: zahl
5652 character(len=5) :: liczba5
5654 integer :: npdb !,ilen
5658 ! real(kind=8) :: var,ene
5662 if (npdb.lt.1000) then
5663 call numstr(npdb,zahl)
5664 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5666 if (npdb.lt.10000) then
5667 write(liczba5,'(i1,i4)') 0,npdb
5669 write(liczba5,'(i5)') npdb
5671 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5673 call pdbout(ee,titelloc1,ipdb)
5676 end subroutine write_pdb
5677 !-----------------------------------------------------------------------------
5679 !-----------------------------------------------------------------------------
5680 subroutine write_thread_summary
5681 ! Thread the sequence through a database of known structures
5682 use control_data, only: refstr
5684 use energy_data, only: n_ene_comp
5686 ! implicit real*8 (a-h,o-z)
5687 ! include 'DIMENSIONS'
5689 use MPI_data !include 'COMMON.INFO'
5692 ! include 'COMMON.CONTROL'
5693 ! include 'COMMON.CHAIN'
5694 ! include 'COMMON.DBASE'
5695 ! include 'COMMON.INTERACT'
5696 ! include 'COMMON.VAR'
5697 ! include 'COMMON.THREAD'
5698 ! include 'COMMON.FFIELD'
5699 ! include 'COMMON.SBRIDGE'
5700 ! include 'COMMON.HEADER'
5701 ! include 'COMMON.NAMES'
5702 ! include 'COMMON.IOUNITS'
5703 ! include 'COMMON.TIME1'
5705 integer,dimension(maxthread) :: ip
5706 real(kind=8),dimension(0:n_ene) :: energia
5708 integer :: i,j,ii,jj,ipj,ik,kk,ist
5709 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5711 write (iout,'(30x,a/)') &
5712 ' *********** Summary threading statistics ************'
5713 write (iout,'(a)') 'Initial energies:'
5714 write (iout,'(a4,2x,a12,14a14,3a8)') &
5715 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5716 'RMSnat','NatCONT','NNCONT','RMS'
5717 ! Energy sort patterns
5722 enet=ener(n_ene-1,ip(i))
5725 if (ener(n_ene-1,ip(j)).lt.enet) then
5727 enet=ener(n_ene-1,ip(j))
5739 ist=nres_base(2,ii)+ipatt(2,i)
5741 energia(i)=ener0(kk,i)
5743 etot=ener0(n_ene_comp+1,i)
5744 rmsnat=ener0(n_ene_comp+2,i)
5745 rms=ener0(n_ene_comp+3,i)
5746 frac=ener0(n_ene_comp+4,i)
5747 frac_nn=ener0(n_ene_comp+5,i)
5750 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5751 i,str_nam(ii),ist+1,&
5752 (energia(print_order(kk)),kk=1,nprint_ene),&
5753 etot,rmsnat,frac,frac_nn,rms
5755 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5756 i,str_nam(ii),ist+1,&
5757 (energia(print_order(kk)),kk=1,nprint_ene),etot
5760 write (iout,'(//a)') 'Final energies:'
5761 write (iout,'(a4,2x,a12,17a14,3a8)') &
5762 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5763 'RMSnat','NatCONT','NNCONT','RMS'
5767 ist=nres_base(2,ii)+ipatt(2,i)
5769 energia(kk)=ener(kk,ik)
5771 etot=ener(n_ene_comp+1,i)
5772 rmsnat=ener(n_ene_comp+2,i)
5773 rms=ener(n_ene_comp+3,i)
5774 frac=ener(n_ene_comp+4,i)
5775 frac_nn=ener(n_ene_comp+5,i)
5776 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5777 i,str_nam(ii),ist+1,&
5778 (energia(print_order(kk)),kk=1,nprint_ene),&
5779 etot,rmsnat,frac,frac_nn,rms
5781 write (iout,'(/a/)') 'IEXAM array:'
5782 write (iout,'(i5)') nexcl
5784 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5786 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5787 'Max. time for threading step ',max_time_for_thread,&
5788 'Average time for threading step: ',ave_time_for_thread
5790 end subroutine write_thread_summary
5791 !-----------------------------------------------------------------------------
5792 subroutine write_stat_thread(ithread,ipattern,ist)
5794 use energy_data, only: n_ene_comp
5796 ! implicit real*8 (a-h,o-z)
5797 ! include "DIMENSIONS"
5798 ! include "COMMON.CONTROL"
5799 ! include "COMMON.IOUNITS"
5800 ! include "COMMON.THREAD"
5801 ! include "COMMON.FFIELD"
5802 ! include "COMMON.DBASE"
5803 ! include "COMMON.NAMES"
5804 real(kind=8),dimension(0:n_ene) :: energia
5806 integer :: ithread,ipattern,ist,i
5807 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5809 #if defined(AIX) || defined(PGI)
5810 open(istat,file=statname,position='append')
5812 open(istat,file=statname,access='append')
5815 energia(i)=ener(i,ithread)
5817 etot=ener(n_ene_comp+1,ithread)
5818 rmsnat=ener(n_ene_comp+2,ithread)
5819 rms=ener(n_ene_comp+3,ithread)
5820 frac=ener(n_ene_comp+4,ithread)
5821 frac_nn=ener(n_ene_comp+5,ithread)
5822 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5823 ithread,str_nam(ipattern),ist+1,&
5824 (energia(print_order(i)),i=1,nprint_ene),&
5825 etot,rmsnat,frac,frac_nn,rms
5828 end subroutine write_stat_thread
5829 !-----------------------------------------------------------------------------
5831 !-----------------------------------------------------------------------------
5832 end module io_config